o 



NASA CONTRACTOR 
REPORT 


NASA CR-129000 


ROCKET ASCENT 6-LIMITED MOMENT- BALANCED 
OPTIMIZATION PROGRAM (RAGMOP) 

By J. T. Lyons, W. S. Woltosz, G. E. Abercrombie, 
and R. G. Gottlieb 

Northrop Services, Inc. 

P. O. Box 1484 
Huntsville, Alabama 35807 


August 1972 
Final Report 


ST 

a 

n z z vc 

P ? £o 

^ Me 
^ ty o 

35 to 
E c 
s n 

r ° 

< c* 


55 
* o 

' n 

rf 

* c 

# T7 


n fv 



c 

*K3 

M O 
to O 
H to 
^ r*j 
> H 
» J 

M Zt- 

O ;s> 

» o 

to 
to a: 

to H - 

c \ 

ce . 

to * ' 

M \ 

to N # 

3 
H 
H 
to 
G 


to 

-nI 

Oj 

I 

Nj 

-a 

CD 


er 


Prepared for 

NASA-GEORGE C. MARSHALL SPACE FLIGHT CENTER 
Marshall Space Flight Center, Alabama 35812 



I Report No 
NASA CR- 129000 

2 Government Accession No 

4. Title and Subtitle 


ROCKET ASCENT G- LIMITED MOMENT-BALANCED 

OPTIMIZATION PROGRAM (RAGMOP) 

7. Author (s) 


J. T. Lyons, W. S. Wottosz, G. E 

. Abercrombie , R. G. Gottlieb 

9. Performing Organization Name and Address 


Northrop Services, Inc. 


P. O. Box 1484 


Huntsville, Alabama 35807 


12. Sponsoring Agency Name and Address 


National 'eronautics and Space Administration r 

Washington, D. C. 20546 



3 Recipient's Catalog No 


5. Report Date 

St 1972 


6. Performing Organization Code 


8. Performing Organization Report No. 

TR-243- 1078 


10. Work Unit No. 


1 1 . Contract or Grant No 
NAS 8-27621 


13. Type of Report and Period Covered 
Contractor 


14. Sponsoring Agency Code 



16. Abttract 


This document describe* the RAGMOP (Rocks t Ascent G- limited Moment- 
balanced Optimization Program) .©puter program for parametric ascent trajectory 
optimization. RAGMOP computes optimum polynomial- form attitude control histories, 
launch azimuth, engine burn-tlmee, end gross liftoff weight for space shuttle 
type vehicles using a eearch-a.cele rated, gradient projection parameter 
optimization technique. The trajectory model available In RAOIOP includes a 
rotating oblate earth model, Cue option of input wind tables, discrete and/or 
continuous throttling for the purposes of Halting the thrust acceleration 
and/or the maximum dynamic pressure, limitation of the structural load indicators 
qe end q6 (the product of dynamic pressure with angle-of-attack and 9idesllp 
angle), and a wide selection of Intermediate and terminal equality constraints. 

Two step-size control schemes in RAGMOP allow the program to rapidly recover 

fron extremely poor nominal "guess" traje lories and Improve the rate of convergence 

over the classical gradient projection me. hod. RAGMOP is designed to run on 

the Unlvac HOT (Exec 8 version) with less than 30K storage required and typical. 

tun times of from two to tan minutes depending upon the quality of the nominal 

(guess) trajectory end the sophistication of the control program desired. 


17. Key Words (Suggested by Author (s)} 

Pa ran** trie Ascent 
Trajectory 
Computer Program 


18. Distribution Statement 

Unclassified - Unlimited 


E. D Geissler, Director 
Aero-Aa trodvnamics Laboratory 



19. Security Qcmif. (of this repnrt) 

20. Security Classif. (of this page) 

21. No. d Pages 


Unclaimed 

Unclassified 

266 


22. P. ice 
NTIS 


WSrC * Pom 3292*1 < April 1972: 






















FOREWORD 


This technical report presents work performed by Northrop Services, Inc., 
Huntsville, Alabama, while under contract to the Aero-Astrodynamics Laboratory 
of the Marshall Space Flight Center (NAS8-27621) . Mr. R. G. Toelle of 
S&E-AERO-GT was the Contracting Officer's Representative for this task. 


#| 






TABLE OF CONTENTS 


Section Title Page 


FOREWORD ii 

% 

LIST OF ILLUSTRATIONS v 

LIST OF TABLES vii 

PART I - ENGINEERING MANUAL 

I INTRODUCTION 1-1 

II GENERAL DESCRIPTION 2-1 

2.1 GENERAL DESCRIPTION 2-1 

2.2 PARAMETER OPTIMIZATION 2-1 

2.3 CONSTRAINTS 2-2 

2.4 PHYSICAL MODEL 2-4 

2.5 OUTPUT OPTIONS 2-5 

III THEORY 3-1 

3.1 DISCUSSION 3-1 

3.2 PHYSICAL MODEL 3-1 

3.3 OPTIMIZATION/RESTORATION 3-36 

3.4 INTEGRATION PACKAGE 3-52 

IV REFERENCES 4-1 

PART II - PROGRAMMER’S MANUAL 

V MAIN PROGRAM FLOW AND OPERATION 5-1 

5.1 INITIALIZATION 5-1 

5.2 TRAJECTORY INTEGRATION AND EVALUATION 5-1 

5.3 CONVERGENCE TEST 5-4 

5.4 PARTIAL DERIVATIVE CALCULATION 5-4 

5.5 PARAMETER UPDATE 5-4 

5.6 PLOT AND TABLE OUTPUT 5-6 

VI SUBROUTINE DESCRIPTIONS 6-1 

6.1 ALPHABETICALLY ORDERED LIST OF SUBROUTINE NAMES 

AND FUNCTIONS 6-1 

6.2 ALPHABETICALLY ORDERED LIST OF SYSTEM SUBROUTINE 

NAMES AND FUNCTIONS 6-2 

6.3 DETAILED DESCRIPTIONS OF SUBROUTINES 6-3 

VII PROGRAM LISTING 7-1 


Ui 



TABLE OF CONTENTS (Concluded) 


Section Title Page 


VIII VARIABLE NAME CROSS REFERENCE OF MAJOR SUBROUTINES 8-1 

PART III - USER’S MANUAL 

IX INPUT /OUTPUT 9-1 

9.1 PROBLEM FORMULATION 9-1 

9.2 ITEMIZED DESCRIPTION OF NAMELIST INPUT 9-7 

9.3 D:- TAILED OUTPUT DESCRIPTION 9-22 

9.4 CONSTRUCTION OF THE DATA DECK 9-53 

9.5 SAMPLE PROBLEM 9-65 

Appendix A - COORDINATE TRANSFORMATIONS A-l 

Appendix B - INTERPOLATION METHODS B-l 

Appendix C - SOLUTION OF MOMENT BALANCE EQUATIONS C-l 

Appendix D - PAYLOAD CALCULATION D-l 

Appendix E - FLYBACK FUEL CALCULATION E-l 

Appendix F - VELOCITY LOSS EQUATIONS F-l 

Appendix G - PAYOFF/CONSTRAINT PARTIAL DERIVATIVES WITH 

RESPECT TO FINAL STATE G-l 

\ 



LIST OF ILLUSTRATIONS 


Figure 


Title 


Page 


2-1 

3-1 

3-2 

3-3 

3-4 

3-5 

3-6 

3-7 

3-8 

3-9 

3-10 

3-11 

3-12 

3-13 

5-1 

9-1 

9-2 

9-3 

9-4 

9-5 

9-6 

9-7 

9-8 

9-9 

9-10 

9-11 

9-12 

9-13 

9-14 


RAGMOP GENERAL FLOW DIAGRAM 2-6 

EQUATORIAL INERTIAL COORDINATE SYSTEM X Y Z 3-3 

LAUNCH PLUMBLINE INERTIAL COORDINATE SYSTEM X Y Z 3-4 

GEOCENTRIC SPHERICAL COORDINATE SYSTEM r 6 3-5 

BODY AXIS COORDINATE SYSTEM 3-6 

ANGLE-OF-ATTACK (a), SIDESLIP ANGLE (6), CHI-PITCH, 

AND CHI-YAW 3-12 

THE EARTH - AN ELLIPSOID 3-18 

TWO-ENGINE EQUIVALENT THRUST COMPONENTS 3-24 

CENTER OF GRAVITY / GIMBAL PLANE COORDINATE SYSTEM 3-25 

COORDINATED TURN MOMENT BALANCE 3-28 

VEHICLE PITCH PLANE 3-31 

FORCES ON THE RELATIVE VELOCITY SYSTEM 3-32 

ANGLE-OF-ATTACK LIMITING WITH COORDINATED TURNS 3-33 

V V x r 3-35 

RAGMOP MAIN PROGRAM ’ JO W DIAGRAM 5-2 

THRUST EVENT PROFILE 9-2 

INPUT BODY AXIS SYSTEM 9-6 

TRAJECTORY BLOCK PRINTOUT 9-26 

PAYLOAD WEIGHT BLOCK 9-27 

TRAJECTORY CONSTRAINT SUMMARY TABLE 9-27 

INFLUENCE COEFFICIENTS TABLE 9-28 

WEIGHTING MATRIX OUTPUT 9-28 

TOTAL I-SY-SY MATRIX OUTPUT 9-28 

DRHO, WDS, GNU OPTIMIZATION PARAMETERS 9-29 

CONTROL STEP VARIABLES PRINTOUT 9-29 

PARAMETER UPDATE AND CONVERGENCE TEST BLOCK 9-30 

QY SEARCH OUTPUT 9-31 

ORBITAL ELEMENT SUMMARY TABLE 9-32 

CHI PARAMETER VALUES 9-33 


v 



LIST OF ILLUSTRATIONS (Concluded) 


Figure Title Page 


9-15 B0PTBL OUTPUT 9-35 

9-16 PAPER OUTPUT FOR B0PPLT 9-47 

9-17 B0PPLT PLOT OUTPUT 9-47 

9-18 SAMPLE CONTROL SETUP FOR MSFC UNIVAC 1108 EXEC VIII SYSTEM . . 9-53 

9-19 SAMPLE DECK SETUP WITH NTABLE * 0 9-63 

9-20 B0PTBL SAMPLE DECK SETUP 9-64 

9-21 B0PPLT SAMPLE DECK SETUP 9-64 

9-22 SAMPLE DECK SETUP USING B0PTBL AND B0PPLT 9-65 

9-23 NAMELIST INPUT 9-66 

9-24 NOMINAL (INITIAL GUESS) TRAJECTORY 9-75 

9-25 NOMINAL TRAJECTORY SUMMARY. BEGIN QY SEARCH 9-79 

9-26 RESULTS OF QY SEARCH STEP. CONVERGENCE TEST PASSED 9-82 

9-27 CONVERGED TRAJECTORY 9-84 



LIST OF TABLES 


Table Title Page 


2-1 CONSTRAINT CODES 2-3 

6-1 FUNCTION LIBRARY FOR CONSTRAINTS 6-4 

9-1 CONSTRAINT CODES 9-19 

9-2 OPTIMIZATION CODING 9-19 

9-3 THRUST EVENT CUTOFF OPTIONS 9-20 

9-4 USE OF LPRINT OPTION 9-21 

9-5 B0PTBL INPUT VARIABLES 9-34 

9-6 ALPHANUMERIC PARAMETER WORDS FOR B0PPLT INPUT 9-45 

9-7 ALPHANUMERIC UNIT WORDS FOR B0PPLT INPUT 9-45 

9-8 NAMELIST CONSTANTS 9-55 

9-9 RAGMOP NAMELIST INPUT SYMBOLS (ORDERED ACCORDING TO USE) - 

TRAJECTORY PHYSICAL MODEL 9-58 

9-10 RAGMOP NAMELIST INPUT SYMBOLS (ORDERED ACCORDING TO USE) - 

PROGRAM CONTROL 9-62 

9-11 RAGMOP NAMELIST INPUT SYMBOLS (ORDERED ACCORDING TO USE) - 

TRAJECTORY AND CONTROL 9-62 



PART I 

ENGINEERING 

MANUAL 



Section I 

INTRODUCTION 


The advent of the space shuttle produced the requirement at NASA, 

Marshall Space Flight Center, for a computer program that could be used as 
an analysis tool to evaluate several different vehicle configurations under 
various launch conditions, A program with the sophistication of the physical 

(D* 

model in the ROBOT computer program was desired, with the additional 

capability of moment-balanced, lifting trajectories in the atmosphere from 

launch to orbital injection. In addition to the moment-balance, several 

additional constraints were required to be met such as: maximum dynamic 

pressure, maximum qa and qB (product of dynamic pressure with angle-of-attack 

and sideslip angle), booster flyback fuel requirements, and the calculation of 

actual payload rather than injected weight. A number of approaches have been 

taken to the ascent trajectory optimization problem, including: calculus of 

variations^, steepest *escent^*^\ min-H^ J \ gradient projection^ 5 ^*^ , 

(9) 

and conjugate gradient methods . The emergence in recent years of accelerated 
gradient projection methods as perhaps the most powerful parameter optimization 
techniques available for highly nonlinear systems^^*^ , ’^*^‘^ led to the 
selection of gradient projection as the optimization scheme in ft/VGMOP (Rccke*. 
Ascent G-limited ^Moment-balanced Optimization Program). The desire was also, 
of course, for a program that was as compact, last, and easy to use as possible. 

The resulting program, RAGMOP, has the capability of computing optimal engine 
burn-times, liftoff weight, launch azimuth, and polynomial-form attitude 
histories including the effects of atmospheric flight from launch to orbit. A 
static moment-balance scheme balances moments using thrust vectoring in all- 
stages of the vehicle. A large variety of constraints, equality and inequality, 
intermediate and terminal, are available as well as a widely variable control 
program. The program occupies less than 30K storage and converged lifting 
moment-balanced trajectory runs have been obtained in less than 3 minutes on a 
Univac 1108 with nominal "guess" trajectories in error as much as 560 meters/ 
sec in velocity, 5 degrees in flight-path angle, and 176 km in radius at injection. 


* Superscript numbers refer to references listed in Section TV. 



Tfifs document Is arranged In three parts plus an appendix. Part I 
includes the first four Sections and may be considered an engineering manua 1 
Part II is comprised of Sections V - VIII and serves the purposes of a 
programmers manual. Part III is made up of Section IX and is sufficient 
within itself as a user’s manual. The Appendices follow Part III. 

Part I - Enginee r ing Manual 


Section 

I 

Introduction 

Section 

II 

General Description 

Section 

III 

Theory 

Section 

IV 

References 


Part II - Programmer’ s Manual 


Section 
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Program Flow and Operation 

Section 

VI 

Subroutine Descriptions 

Section 

VII 

Program Listing 

Section 

VIII 

Variable Name Cross Reference 



Part III - User’s Manual 

Section 

IX 

Input/Output 
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Section II 

GENERAL DESCRIPTION 

2.1 GENERAL DESCRI PTION 

The RAGMOP (Rocket Ascent G-limited Moment-Balanced Optimization Program) 
computer program calculates the optimal value** of a set of paramecers which af 
feet multistage rocket ascent trajectories. rhe parameter optimization is per- 
formed using a search-accelerated gradient projection technique, which includes 
the capability of satisfying a large number of intermediate and terminal con- 
straints (end conditions) . The flexibility and speed of the progr^ , combined 
with a highly sophisticated physical model, make it 2 desirable tool for the 
analysis and design of space shuttle and other ascent rocket vehicles. ..a 
addition, output options include tables and plots which are suitable for reports. 

2.2 PARAMETER OPTIMIZATION 

The parameter optimization scheme used by RAGMOP is a search-accelerated 
gradient projection method. This method requires complete parameterization of 
all control variables , and the solution obtained is therefore optimum only to 
the degree attainable with the chosen parametric form. 

The control parameters in RAGMOP describe: (1) the lift-off weight of the 
vehicle, (2) the duration of a number of engine burn times (thrust events), 

(3) the launch azimuth, and (4) the vehicle pitch and yaw attitude histories. 

The vehicle chi-pitch (x p ) and chi-yaw (x^) attitude histories are gi> 
the form of polynomials in time, l*e.: 

Xp " Xp + a x (t-t o ) + a 2 (t-t o ) 2 + • . • + a n (t-t o ) n 

0 2 
and x Y * X Y + b i ( b_t 0 ) + b 2^ t_t 0 ^ 
o 

Separate polynomials are used for each stage, with currently up to a fourth- 
order polynomial available for the first stage Xp» and second-order polynomials 
for the second stage Xp and the x Y program of both stages. The Xp program may 
be continuous or discont'nuouc: at staging, while the x Y program is always con- 
tinuous. The validity of the polynomial form for the Xp and x Y attitudes angles 
(See Section 111 for mere information concerning the attitude control) is evi- 
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denced by the close agreement between polynomial forms and solutions obtained 
using variational methods. 

In addition to the optimized parametric form x mentioned above, RAGMOP 

P 

also allows the use of an angle-of-attack profile for the pitch attitude con- 
trol of a portion of the first stage f ght. By the use of an input flag, the 
user may specify three angle-of-attack control options: (1) zero aerodynamic 

normal force, (2) zero angle-of-attack, or (3) angle-of-attack as a function 
of Mach number. The angle of attack control will be used from the end of a 
tilt-over maneuver to staging. The tilt-over will be performed between the 
end of the lift-off (vertical rise) phase and the beginning of angle-of-attack 
control at some time specified by the user. The tilt-over consists of an op- 
timized polynomial Xp control program of the same form as in the complete stage 
when angle-of-attack control is not used* 

Revision I adds the capability of enforcing coordinated turns during the 
first sta^e. Either positive or negative angle-of-attack is used for this 
option, the algebraic sign determined separately for each thrust event. This 
option may be useu simultaneously with the angle-of-attack options mentioned 
above. 

2.3 CONSTRAINTS 

RAGMOP is extremely flexible in terms of the constraints flowed on the 
trajectory. 

Table 2-1 lists the equality constraints available in the program, which 
may be enforced both at orbital injection and at staging. In addition to the 
equality constraint c£ Table 2-1, the relative velocity at staging may be used 
as a cutoff criteria for the last thrust event oi the first stage. This cut- 
off criteria will be satisfied regardless of che r.arameter values (if possible) 
and does not enter into the parameter update equations. 

Als av 'liable in RAGMOP are several inequality constraints, namely: (1) 

the product dynamic pressure with angle-of-attack (qu, , ( 2 ; the product of 
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tlvu-imit pressure with sideslip angle (qB), and (3) the maximum thrust accelera- 
tion ol the vehicle (g-limit). 


The qa and qB constr ' its are enforced by reducing a ar.d/or B to produce 
the maximum acceptable values whenever they are exceeded. This results in 
temporarily overriding the Xp and/or Xy polynomials until such time as the 
q i and qB produced by the polynomials is acceptable. 

Table 2-1. CONSTRAINT COOES 

The codes contained in this table are input into KCDPHI fid KCDRES to 
designate the payoff and the intermediate and terminal constraints desired 
for the trajectory. The appropriate values desired for these constraints must 
then be Input into PS1REQ and PSIRST . 


CODE NUMBER 

UNITS 

CONSTRAINT OR PAYOFF 




1 

KG 

Payload 

2 

M/SEC 

Inertial velocity 

3 

DEG. 

Inertial flight path angle 

4 


Radius 

5 

m 2 /sec 

Energy 

6 

M /SEC 

Angular momentum 

7 

DEG. 

Inertial longitude 

8 

DEG. 

Inertial heading angle (+ East from South) 

9 

DEG. 

Colatitude 

10 

DEG. 

Inclination 

11 

DEG. 

Line of modes 

12 

M 

Seml-latus rectum 

13 


Eccentricity 

14 

SEC 

Total bum time 

15 

LB /FT 2 

Maximum dynamic pressure 

16 

DEG. 

True anomaly 

17 

DEG. 

Argument of perigee 

18 


Reserved for future use 

19 


Reserved for future use 

20 

NM 

Flyback range 


The acceleration (g) limiting can be enforced in two ways: (1) continuous 

throttling may be employed to "ride” the g-limit, if the actual vehicle engines 
allow this, or (2) discrete throttling (shutting down one or more engines or 
reducing the output of all engines in discrete amount! ) may be used. Continu- 
ous throttling is performed whenever the acceleration reaches the desired limit- 
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ing value, at which time a gradual continuous reduction of thrust is accomp 1 ished 
until the er * of the thrust event. Discrete throttling is accomplished bv ini- 
tiating a new thrust event with a lower fixed (except for the exit plane pres- 
sure difference) thrust whenever the desired acceleration limit is reached. 

The following thrust event will begin with an acceleration less than the limited 
value, with the acceleration increasing with time until its limit or another 
thrust event cutoff criteria is reached to begin the next thrust event (see 
Section IX for a more complete explanation of thrust events). 

2.4 PHYSICAL MODEL 

The physical model employed in RAGMDP has been designed to be as realis- 
tic as possible and yet allow the program to remain within the attempted goals 
of 32K storage and 3 minute run time. The result is a highly sophisticated 
trajectory model which provides the user with a number of options in setting 
up rhe trajectory run. 

The geophysical model presents the options of a spherical or oblate, ro- 
tating or nonrotating earth with a tabulated Patrick Reference Atmosphere (1963 
version) and the ability to specify input wind directions and speeds as alti- 
tude l unctions. This last option allows the program to bias the trajectory 
profile to either take advantage of, or to minimize the losses from, winds at 
various altitudes. 

RACMOP also provides the user with the option of thrust vectoring to pro- 
duce zero total moment on the vehicle. A two engine equivalent thrust model 
centered around the actual vehicle thrust centroid is used to reduce the com- 
putation time required for this option. Thrust components are computed to bal- 
ance aerodynamic moments using small angle approximations resulting in negli- 
gible error for small gimbal angles. The error, such as it is, will be in the 
form of slightly unbalanced moments (less than about 1-1/2 percent error for 
gimbal angles of 10 degrees). The total thrust aval able remains intact. Re- 
fer to Section III and Appendix C for a complete discussion cf the moment bal- 
ance scheme. 
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2.5 OUTPUT OPTIONS 


RAGMOP includes two special output subroutines which provide the user a 
set of tables and/or plots summarizing the converged trajectory (solution). 

The tables * which use a fixed format, and the plots, are suitable for publi- 
cation. The output plots are produced on the CALCOMP plotter, and therefore 
use of this option is restricted to systems which iiave the CALCOMP plotter 
available. Plots are produced for any variable versus any other, in any units 
the user desires. 

Figure 2-1 presents a macro-flow diagram of the RAGMOP computer program. 
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RACHOP CHIRAL HOW DIAGRAM 



Figure 2-1. RAGMOP GENERAL FLOW DIAGRAM 
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Section III 
THEORY 


DISCUSSION 

The RAGMOP computer program is designed to solve the rocket ascent 
optimization problem including the effects of trajectory shaping, engine burn 
times, liftoff weight, and launch azimuth* The optimization is performed 
subject to a number of constraints both during and at the end of the trajectory. 

a order to perform this function, three main requirements have been satis- 
fied by the computer program; (1) a physical model has been programmed into 
the equations of motion which represents the actual flight of the vehicle as 
closely as possible, (2) the optimization method attempts to be rapid in 
terms of computer time without unnecessarily restricting the physical model, 
and (3) the scheme used to integrate the equations of motion is rapid and 
flexible with a minimum amount of error in the integration. This section 
presents a detailed description of the three areas just mentioned. 

3.2 PHYSICAL MODEL 

The RAGMOP computer program has been designed to include as sophisticated 
a physical model as possible subject to the computer run time and storage goals 
of three minutes and 32k, respectively. A three-dimensional trajectory model 
with a static moment balance is used with tabulated atmospheric, aerodynamic, 
and center-of-gravity data. The atmosphere model includes a spline-interpo- 
lated Patrick Reference Atmosphere (1963) and input wind tables which allow 
wind directlu and speed to be specified at up to 25 altitudes. Aerodynamic 
data is .'so spline- interpolated* and is input for both stages in the form of 
force and moment coefficients and their angle-of-attack or sideslip angle 
di ivatives, at up to 25 Mach numbers. The earth model includes a rotating 
atmosphere and an oblate (Fischer ellipsoid) gravitational model. A complete 
description of the force and moment equations and the computation of their 
component parts is presented In the following paragraphs. 


* See Appendix B for further information concerning interpolation methods. 



3.2.1 Generalized Equations of Motion 

The equations of mot .on in any inertial reference frame are: 

X(t) = F(t)/m(t) 
t 

• Vo + J * 

0 

t t 

*(‘) ■ Vo + f * dt + , 

0 0 

where 

X is a three-dimensional position vector, 

F is the total force acting on the vehicle, 
m is the instantaneous mass of the vehicle, 

and 

t is time measured from some reference time t«0. 

The moment equation used for the static (3D) moment balance is: 

M « m a + Mr ® ° 

where 

M is the total moment acting on the vehicle, 

M a is the total aerodynamic moment, 

A 

and 

Hp is the total thrust moment. 

3*2.2 Coordinate Systems 

Several coordinate systems are used in the RAGMOP computer program, and 
an understanding of these systems (and the transformations which allow 


t 

| X dt dt 

0 
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changing from one system to another) is essential to a thorough comprehension 
of the equations of motion. The five coordinate systems used in RAGMOP are; 
(1) the equatorial inertial system (2) the launch plumbline inertial system, 
(3) the spherical geocentric system, (4) the body axis system, and (5) the 
relative velocity system. 

3. 2. 2.1 Equatorial Inertial System . The basic reference coordinate system 
in RAGMOP is the equatorial inertial geocentric cartesian coordinate system 
shown in Fig. 3-1. This coordinate system has the Y axis pointing north, the 
X and Z-axes in the equatorial plane, and the Z axis contained in the longi- 
tudinal plane of the launch site. 


NORTH POLE 

St 

y 



*(GYRO reference release ) time at which guidance system ie released 

Figure 3-1 . EQUATORIAL INERTIAL COORDINATE SYSTEM X Y Z. 


«* 

5 

t 

* 

I 



i 
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3. 2. 2. 2 Launch Inertial Plumbline Coordinate System . The launch inertial 
plumbline coordinate system is the system from which the Xp and Xy attitude 
angles are defined and in which the equations of motion are written. The 
y axis of this system is parallel to the launch site gravity vector (plumb- 
line) but, for an oblate earth, does not pass directly through the launch 
site. The x axis is pointed in the direction of the launch azimuth, and the 
z axis forms a right-hand system. The origin of the xyz system is at the 
center of the earth. 



Figure 3-2. LAUNCH PLUMBLINE INERTIAL COORDINATE SYSTEM XYZ 
3.2.2. 3 Spherical Geocentric Polar Coordinate System . The spherical 

A A A 

geocentric polar coordinate system moves with the vehicle. The $ r 6 axes 
of this geocentric system point in the directions of increasing <t>, r and 9, 
respectively, where <J> is measured from the equatorial inertial Z axis to the 
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plane containing r and the equatorial inertial Y axis, and 0 is measured from 
the equatorial inertial Y axis to the r vector. (See Figure 3-3). 


fi 

I vwv, 


m 



Figure 3-3. GEOCENTRIC SPHERICAL COORDINATE SYSTEM «re 

3. 2. 2. 4. Body Axis Coordinate System. The body axis coordinate system, shown 
in Figure 3-4, is the system in which the aerodynamic and thrust forces and 
moments are calculated. The forces are then transformed into the launch 
inertial plumbline system to determine accelerations for integration in the 
equations of motion. The body axis system is defined with the X' ,Y', and Z* 
axes such that the vehicle longitudinal axis is parallel to the Y 1 axis, the 
X 1 axis is positive "downward" with respect to the vehicle, and the Z'axis 
points in the direction of the right wing. This is such that, on the launch 
pad, the vehicle X', Y f , and Z 1 axis are parallel to the launch inertial 







Note that for input purposes only, the body axis system is as shown in 
the insert of Figure 3*4. Input data uses the input body coordinate system. 

The equations of motion use the x'y'z* system described above. (Data is 
converted in the input subroutine). 

3.2. 2.5 Relative Velocity Coordinate System . The relative velocity coordinate 
system is used whenever coordinated turns are specified for the RAGMGP trajec- 
tory (see paragraph 3. 2. 6. 2). This nonorthogonal coordinate system allows the 
use of simple relationships for the aerodynamic and thrust forces on the vehicle, 
resulting in a savings of computer time when the coordinated turn option is used. 
The coordinate system is defined by the plane of the relative velocity vector 
(V ) and the vehicle longitudinal axis (y f ). Vectors in the V_ and y f directions, 

K R 

which are not generally orthogonal, are used to specify the forces in the plane, 

which are the only aerodynamic or thrust torces on the vehicle when using co- 
ordinated turns. When this system is used, the unit vectors of the body axis 
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system are computed in terms of the inertial plumbline unit vectors 9 from which 
the inertial attitude angles Xp» Xy* and are found* 


3* 2. 2. 6 Transformations * A vector in any coordinate system may be transformed 
into any other coordinate system by premultiplication with the proper trans- 
formation matrix* This transformation can be performed as: 

x ® Ax where 

x' is a vector in the x'y'z* coordinate system, 
x is a vector in the x,y,z, coordinate system, 
and A is the matrix which transforms x into x 1 , i*e. the matrix that computes 
the components of x in x y z f so that x can be rewritten as x (x and x are 
the same vector since transformation does not change the vector, but only the 
coordinate system in which it is written).* 

Equatorial Inertial to Launch Plumb line Inertial . The launch plumbline 
inertial system is obtained from the equatorial inertial system by first 
rotating about the equatorial inertial X axis an angle, 0 = 90° - lat 

launch site* 

and then about the y axis (formed from that rotation) an angle, <|i = -(AZ, 

launch 

90°). The transformation matrix from the equatorial inertial XTZ system to 
the launch plumbline inertial xyz system is, then*: 


A Xx 


cos(90-Az) 

0 

sin(90-Az) 


0 

1 

0 


-sin(90-Az) 


co8(90-Az) 


1 0 

0 cos0 

0 -sin0 


0 

sin0 

COS0 


(noting that sin(90-lat) = cos(lat), cos(90-lat) 
cos Az, and cos(90-Az) = sin Az). 


sin(lat) , 8in(90-Az) 


'■~Tt - ’V : 





sin Az 
0 

-cos Az 


sin0 cos Az 

COS0 

sinQ sin Az 


-cos Az 
-sin0 

cose sin Az 


*See Appendix A for further information on coordinate transformations. 
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where lat « geodetic latitude of the launch site and Az » launch azimuth. 


This transformation is used for computing the earth rotational and 

gravitational acceleration components in the launch plumbline system from 

the equatorial system. The reverse transformation uses the transpose of the 

T 

above matrix, A xX - [A^] * 

Body Axis to Launch Plumbline Inertial System. The attitude angles chi- 
pitch ( Xp ) and chi-yaw (x^) describe the orientation of the vehicle with 
respect to the launch plumbline system. The vehicle attitude is obtained by 
first rotating about the launch plumbline z axis the angle -x » and then 
rotating about the vehicle x 1 axis the angle x^* Since, thrust and aerodynamic 
forces on the vehicle are determined in the body axis system, but their 
resultant accelerations (added to the gravitational acceleration) are integrated 
in the launch plumbline system, the transformation matrix from the body axis 
system to the launch plumbline system is required. This is found by noting 
that the transformation involves first rotating about the body x f axis an 
angle and then about the launch plumbline z axis an angle x p so that, 
from the body axis x f y'z f system to the launch plumbline xyz system we have: 


cos X 


sin * 


Xp 0 


Noting that cos - 


cos(-Xy) 

-sin(-Xy) 


cos x_ and sin -) 


sin(-x y ) 

cos(-x y ) 


cos x p 
-sin x p 


sin x p cos Xy 
cos x p cos Xy 
sin Xy 


-sin x p sin Xy 
-cos x p sin Xy 
cos Xy 


Equatorial Inertial to Spherical Geocentric Polar Coordinate System . The 
use of the spherical geocentric polar coordinate system makes the determination 
of a number of parameters more convenient, hence the transformation matrices 


/ 




* 


9 


from the t. y stem to the equatorial and vice versa are necessary. The spherical 
geocentric polar system is obtained by rotating first around the Y axis an 
angle 4> and then about the <t> axis and angle 0, so that: 


1 

0 

0 


cos<(> 

0 

-sin(j) 


COfc i> 

0 

-sin£ 

0 

COS0 

sin0 


0 

1 

0 

= 

sin^siiiO 

cos0 

cosisin6 

0 

-sin0 

COS0 


sin$ 

0 

COS<|> 

— * 


sinj)co30 

-sinb 

COS$COS0 


Note, however, (see Figure 3-3) that: 

Y X 

cos0 = — , tan <J> = — , 
r Z 


JTTS = r sin0, 


so that sin 


>1 x 2 + z 2 ~ 


s inti = 


r sin8 



and 


cos<J> 


r sin0 




so that 




\j x 2 + 


X 

r 


XY 




2 2 
+ Z 


Y 

r 


i 


2 2 
X + Z 


\J X 2 ■¥ Z 2 


Z 

r 


r J7~¥~P 


The reverse transformation, from to X is accomplished using* A 


<f>X 


t A X<^ 


Launch Plumbline to Spherical Geocentric Polar Coordinates > The trans- 
formation from the launch plumbline system to the spherical geocentric polar 
coordinate system can be obtained by first transforming to the equatorial 
Inertial system and then to the spherical system, so that 


Correction page for page 3-9 



A X* “ A x<t> A xX . 

The above matrix product yields (letting A denote Ay and D denote A ) : 


D 11 

S 

(A 22 Z-A 32 Y)/rsine 

°21 

= 

(A 32 X-A 12 Z)/rsin6 

°31 

s 

(A 12 Y " A 22 X)/r8in0 

X 

n 

33 


U 12 


r 


= 

Y 

22 


r 


s 

Z 

32 


r 

°13 

= 

(D 12 cos0-A 12 )/sin0 

°23 

as 

(D 22 cos0-A 22 )/sin0 

°33 

S3 

( D 32 cos 0-A 3 2)/sine 


The transformation from <{>r 9 to xyz (plumbline) is accomplished using the 

T 

transpose of the above matrix; A ^ . 

2.2.3 Forces and Accelerations 

The forces acting on the vehicle in the RAGMDP simulation are the 
thrust, aerodynamic, ar. gravitational forces. The aei dynamic and thrust 
forces are computed in the body axis system and then transformed to the 
launch plumbline system wherein the accelerations are added to the gravitational 
acceleration components for integration in the equations of 

3.2.3. 1 Aerodynamic Forces. The aerodynamic forces acting on the vehicle are 

the axial, normal, and side forces. These forces are found by determining 

the proper axial, normal, and side force coefficients for the vehicle at the 

current Mach number, and then mi *tiplying these coefficients by the product 

2 

of the dynamic pressure (Q ® 1/2 p V ) and the aerodynamic reference area 
(S) , such ' * “ 


[ ] 


T 


denotes matrix transponee . 
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* 

I 


V 


f aero * c axial* q * s> 

AXIAL 


'““Wl ' C,l0 “* u - 


*Q*s» 


^AERO * S»IDE*^*^' ax * a * f° rce is positive aft (-y*), t e 

SIDE 

normal force positive upward on the vehicle (-x') and the side force positive 
toward the right wing (+ 2 '). The aerodynamic coefficients are determined from 
the Mach number and angle— *of '—attack. or sideslip angle by means of a spline 
interpolation* in the input aerodynamic coefficient tables. The dependence of 
coefficients on angle— of —attack (a) or sideslip angle (8) is assumed linear 
at each Mach number, thus: 


C AXIAL = °A + (C A } (a) 
o 01 

SoRMAL * S + (Cg > 

o a 


c side " (e) 


The axial aerodynamic c orce is also affected by the base drag term, which is 
detet^ned using a linear interpolation table lookup from the input tables of 
base axial force and altitude. Thus 


W ■ F <«-TITUDE) 


and then F 


AERO,,,,.., “ C AXIAL ^ ‘ F BASE (The ne « ative si 8 n is used since the 
AXIAL 


axial force (drag) and the base pressure force are opposite in direction. 


The F 


AERO 


AXIAL 


since F 


AERO. 


will be subtracted from the y* thrust component 
is positive in the -y ' direction) . 


AXIAL 

Thus, the total aerodynamic forces in the body axis x'y'z' system are: 

F a = -FAN - - (Cjj + a)QS 

x’ o a 

„ . -EAA . - (C ♦ 0 WO. f BASE 

y o a 


* See Appendix B fc • further in forme v'.on 


.ticm methods . 


1 


i 


? 

2 





'A z , * F SIDE “ (C y s 6)QS 


These forces are transformed Into the launch pluabline system using the 
transformation matrix A , of Paragraph 3. 2. 2. 6, so chat 

X X 


v> 


[vj 


fv' 

V 


The Mach number used in the spline interpolation lookup of the aerodynamic 
coefficients is found by: 

M * VR/a 

where VR = relative velocity (velocity of vehicle with respect to the 
atmosphere) 

and a = speed of sound at current altitude from the spline-interpolated 
atmosphere routine PRA63 

The relationship between the relative velocity vector and the body axis 
system determines the angle-of-attack and sideslip angle of the vehicle 
(See Figure 3-5). 



Figure 3-5. ANGLE-OF-ATTACK (a), SIDESLIP ANGLE ( 0 ), CHI-PITCH, AND CHI-YAW 
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The relative velocity VR is computed as the difference between the plumb- 
1 ine inertial velocity and the transformed velocity required by an object in 
order to remain over a given point on the rotating earth (at the same altitude 
as the vehicle' plus the wind at the current altitude. This results in 


u > -D ) 0 


rft sinVK ~ a 32 y ^ n e 


u -<« 32 « - 
y-(a 12 y - a 22 x » e 


V I cos A E + sin A N\ 

w V z w z w / 


where 


(w,u,v) are the relative velocity components in the plumb line system, 

(w,u,v) are the actual inertial veloci y components of the vehicle in 
the plumbline system. 


( a 22 z-a 32^ ^ a 32 X ~ a 12 Z ^ ^ a 12 y ~ a 22 X ^ ^e are the inertial 

velocity components in the plumbline system required to remain stationary 
over the subpoint on the rotating earth. 

D - A^ x (see Paragraph 3. 2. 2. 6), 

A, is the wind azimuth (+ CV/ from north) 

z w 


V is the wind speed 

w 


E is a unit vector in the east direction 


U x r 
Ik x r| 


N is a unit vector in the north direction (r x E) 
The relative velocity VR is then equal to: 


VR =\/j^ + u^ + v^ 


The dynamic pressure Q is obtained from the relation: 


where 


Q - 1/2 p (VR) 4 

% 

o * atmospheric density from the spline- interpolated Patrick 1963 atmosphere 
routine . 


VR * relative velocity described above* 


* ***1 





N 


3.2. 3.2 Thrust Forces , Thrust forces in RAGMOP are computed in two ways; 
(1) the SRM thrust is based on input values of sea level thrust corrected 
Lor the atmospheric pressure at the current altitude, (2) the liquid engine 
thrust is based on vacuum thrust corrected for current atmospheric pressure. 


The thrust calculation for the SRM*s is as follows: 


T = (T + A^P.. - P AM ))N., Kir 
SL L SL AM bNG 


where 




and 

N 


ENG 


total thrust available, 

sea level thrust per SRM engine from spline-interpolated tables 

exit area per engine, 

sea level static pressure, 

ambient pressure at current altitude, 

number of SRM engines used for the current thrust event. 


The liquid engine thrust calculation utilizes input values lor vacuum 
thrust and corrects for the exit plane pressure differential: 


T = (T 


VAC 


a e p am )n eng 


where 

T ^ * vacuum thrust per engine 
v A v 

Acceleration-Limited Thrust . Two methods exist in RAGMOP for limiting 
the acceleration of the vehicle (g-limit). One method available is to simply 
create a new thrust event any time the acceleration reaches the desired lim- 
iting value. The new thrust event can have fewer engines* or all engines 
operating at a reduced thrust level. The second method utilizes continuous 
throttling to hold the acceleration at the limited value. The approximation 
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is made that the acceleration of the vehicle due to the thrust force above will 
be equal to the total acceleration. Thus, in the x'y'z' body axis coordinate 
system. 


T 


8 limit 

m 


if g > 


8 limit 


where * maximum acceptable longitudinal acceleration, 

m * instantaneous mass of the vehicle, and 
T ■ total vehicle thrust. 


Ja order to enforce the continuous throttling g-limit constraint, the value of 
the thrust acceleration is checked by the integration package until the limit 
is reacVed. The integration package isolates on the exact time that g-limit 
occurs* and then either performs thrust event initiation (discrete or step 
throttling) or begins continuous throttling as described above. Note that for 


the SRM booster, only the main (liquid) engines are throttled unless all 
engines are used in the moment balance (SRM*s and main engines). 

3. 2. 3. 3 Accelerations . The total acceleration components in the 4 launch plumb- 
line inertial coordinate system are simply 


x 

y 

z 


F 

— + g 

m & x 

F 

m s y 

F 

— + g 

m 


where F , F , F are the total forces on the vehicle due to thrust and aero- 
x* y z 

dynamics and g x , g^, are the gravitational accelerations in the launch 
plumbline coordinate system (see Paragraph 3.2.4). Expanding the above equa- 
tions we obtain (where A • comes from Paragraph 3.2. 2.6) 

x x 


'x' 




f T . 

X 

akro normal 



y 

* - 

A , 

X X 


T * 

y 

■ r “«w 

> + l 

00 

k* - 


m. 


k T • 
z 

+ f side > 




*See Paragraph 3.4 for more information concerning the integration package. 





Note that and F 

^NORMAL 

directions, respectively. 


AXIAL 


are positive in the negative x* and y* 


3.2.4 Moss Calculation 

The instantaneous mass of the vehicle must be known in order to calculate 
the accelerations in the inertial reference frame. Fuel expended and weight 
jettisoned must be known in order to perform the mass calculation. Fuel for the 
liquid engines is integrated whereas SRM fuel overboard is determined using a 
spline- Interpolated table Lookup. The Jettison weLghts are discrete weight 
drops that can occur only at the end oi any thrust event. This is necessary 
since a weight jettison produces discontinuities in the state derivatives and 
the integration scheme must "restart" when variable step methods are used. 


The mass at any instant * XM» is given by 


XM = XM1AD + XMAUG 


where rt 

XMIAD ® XMI + J -XMDOT 
t=0 

0 

XMl t=() * WZERO - XMAUG t-() 


dt 


and 


NVNT 

XMAUG. , * l WJET (1) 
t '° 1-1 

WJET (I) = jettison weight per thrust event 

NVNT = total number of thrust events 

ITHR-1 

XMAUG = XMAUG - E WJET (I) 

t=0 I=1 

ITHR = number of current thrust event 


When the SRM booster option is used, SRM mass overboard is found from a 
spline interpolation of time-dependent values. The effect of this weight loss 
is included in the mass calculation by using the equation 
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XMIAD = XMI + SRP 

where 

SRP * SRMPRP - SRMDWT 
SRMPRP = total SRM propellant 
SRMDWT = SRM mass overboard 
and 

[' 

XMI * XMI. n + I - XMDflT dt. 
t=0 * 

0 

Thus as before, 

XM = XMIAD + XMAUG. 

When the last thrust event using SRM's is completed, the empty cases are jet- 
tisoned as one of the WJET(I) and the total fuel expended, SRMPRP, is sub- 
tracted from XMAUG, so that the mass equation remains 

XM = XMIAD + XMAUG 


3.2*5 Geophysical Model 

RAGMOP uses an oblate earth model (Fischer ellipsoid) with a preset flat- 
tening coefficient of f = 1/298.3 (which may be changed by input). (See 
Figure 3-6). 

The radius of the earth using this model at any given geocentric co- 
latitude 9 is given by: 


R(9) = (l-f)R / 
e 




2 2 2 
sin 9 + cos 9 


The derivative of radius with respect to colatitude 0 (required in order 
to calculate the time at which maximum dynamic pressure occurs) is: 


dR(9) R (9) f (2-f)sin8cos9 
d9 R 2 (l-f) 2 
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Figure 3-6. THE EARTH - AN ELLIPSOID 
3.2. 5,1 Atmospheric t roperties . The atmospheric model in RAOMOP is a 

( 15 )* 

Patrick Reference Atmosphere, 1963 version , using a spline-interpolation 
technique to obtain the properties of pressure, density, and speed-of-sound, 
and their altitude derivatives* This atmosphere rotates with the earth at 
an angular velocity il = 7.2921158 x 10 ^ rad/sec, resulting in no wind over 
the surface of the earth* Winds may be added as functions of altitude so 
that the trajectory shaping performed by the program is biased to include the 
wind effects. 

The initial inertial velocity of the vehicle (at liftoff) will be the 
velocity produced by the earth’s rotation. The relative velocity at time = 0 
will be zero unless input wind tables specify some ground wind direction and 
speed. 


3.2. 5. 2 Gravitational Model . ^ The Fischer ellipsoid gravitational model 
is used unless the user specifies a spherical earth model in the input. The 
gravitational potential function used for the oblate model is given by: 


* Superscript numbers refer to references in Section IV. 




\ 


u R r 

U(r,0) * 7 s U + Y (-^) 2 ( 1-3 cos 2 0) + | (-^) 3 (3-5 cos 2 0)cosO 


+ ^ (~)^ (3-30 cos 2 0 + 35 cos^0)] where 


CJ » 1.62345 x 10'' 
H = -0.575 x 10‘ 5 
DJ « 0.7875 x 10' 5 


R g ■ earth equatorial radius * 6378165 m 

u * Product of universal gravity constant and earth mass « 3.986032 x 

6 i/»14 3. 2 

10 m /sec . 


Each of the above may be changed by input if desired by the user. 

The components of the gravitational acceleration vector in the launch 
plumbline inertial system are calculated as the first partial derivatives 
of the potential function with respect to the plumbline coordinate axes. Thus, 


s-S 


3U I dr I 3U ( 30 

3r | 3yf 30 3y 


These equations may be rearranged into the form: 


G 11 { y f" G To { a 22 


where G, 


V R 7 , R , 7 

-f [ 1+CJ (*— ) * ( l-5cos o) + H(-V(3-7co8^0)cose 
r r r 


+ DJ(-^) A (-| - (6-9 cos 2 0)coa 2 0)l 
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v. 


G To " ~f [2CJ(~)“ cost! - H(-jr) J (j - 3 cos'o) 


R 


e \ 3 ,3 


R 


+ DJ (~~) ^ (~ y - 4 cos^0) cos0] 


and 

is the second column of the transformation 
matrix from the launch plumbline system to 
the equatorial inertial system (see Paragraph 3. 2. 2. 6). 

This general form is the same as that used in the Saturn V flight computer. 



In the event a spherical earth is specified the above relations reduce to: 



So- 0 - 






3.2,6 Moments 

RAGMOP employs a static (3-D) moment balance scheme using thrust vectoring 
to balance aerodynamic moments In order to more accurately model the actual 
vehicle performance. Solutions obtained without moment -balancing assume the 
entire thrust to be directed through the vehicle center-of-gravity, and are 
suitable for symmetric vehicles such as the Saturn V. However, for a non- 
symmetric vehicle with offset center-of-gravity and thrust centroid (as in 
space shuttle configurations) the results obtained without moment balancing 
can differ greatly from those obtained with the moment balance enforced. 

Moment balancing requires redirecting the vehicle thrust vector so as to 
produce no net moments acting on the vehicle. Thus, the thrust component 
acting through the center-of-gravity (and hence the total vehicle acceleration) 
will be less than that obtained for the nonmoment-balanced case* The moment 
balance scheme employed in RAGMOP makes use of two approximations: (1) the 

controllable engines are lumped into a two-engine equivalent thrust model, and 
(2) a small angle approximation is used in the solution for the required gimbai 
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angles. These approximations, which yield a slight residual unbalanced moment, 
serve to reduce the run time of the program considerably. 




3.2.6. 1 Aerodynamic Moments . The aerodynamic moments are found in a 

manner similar to the aerodynamic forces. Moment coefficients are determined 
using spline interpolation for the tabulated values of the moment coefficients 
and their angle-of-at tack/sideslip angle derivatives as functions of Mach 
number. The coefficients used are: 

CMO - pitching moment coefficient at zero angle-of-at tack, positive 
about z’ according to the right-hand rule 

CMALP - partial derivative of pitching moment coefficient with respect 
to angle-of-at tack. 

CLBETA - partial derivative of rolling moment coefficient with respect 
to sideslip angle. Rolling moment is positive about y f 
according to the right-hand rule. 

CNBETA - partial derivative of yawing moment coefficient with respect 

to sideslip angle. Yawing moment is positive about x' according 
to the right-hand rule. 


The total aerodynamic moment about the center* of-gravity is obtained 
by determining the total moment coefficient about the aerodynamic reference 
point for which the coefficients are specified (and which is not, in general, 
coincident with the center-of-gravity) and then translating these moments 
to the center-of-gravity. The aerodynamic moments about the body axes at 
the center-of-gravity are then determined by multiplying the coefficients 
by the product of dynamic pressure, aerodynamic reference area, and aerodynamic 
reference length. Thus, 

°M , * C M + C M ° 
z o a 





about the aerodynamic reference point. 


and 
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V 

eg 




'eg 


= °M ‘ C N (y cg " y REF^ ^ 11 REF + C Ax* X cg ~ X REF )/£ REF 
z p 

" C M , + C S <X cg “ X REF )/a REF 

y 


°M , * Si , " C S (y cg " y REF )/J, REF 


eg 

where , C^, and are the normal, axial, and side force coefficients of 
P 

Paragraph 3. 2. 3.1, and £ R g F is the aerodynamic reference length of the vehicle. 


The actual aerodynamic moments about the eg are then found by: 


\ , ’ Si , <i' 2 °™ 2 >s Vf ' Si , 0 s ‘ref 

X X X 

eg eg eg 


M i. , " Si , ^ £ ref 

y cg y cg 

V =C M, QS *REF 

z' z 

eg eg 

3. 2. 6. 2 Thrust Moments . The moment balance scheme employed in RAGMOP computes 
the thrust components required to balance moments using a two engine equivalent 
thrust model for the controllable engines. The moment balance computations 
are carried out in ore of several different ways, depending upon the option 
desired by the user. Three options for SRM engines are available: (1) SRM 

engines thrusting through a fixed point in the vehicle x' - y' plane, (2) 

SRM's tracking the c.g., and (3) SRM' s balancing moments along with the main 
(liquid) engines. In addition, the coordinated turn option may or may not be 
used. With coordinated turns, only pitch plane moment balancing is required. 


When SRM's are not used for moment balancing, positions of the two equi- 
valent engines are determined by summing the total liquid engine gimbal point 
X' and Y' locations and dividing by the total number of engines to obtain the 
average X' and Y' gimbal points, and by summing ail of the positive Z' positions 
and dividing by half the number of engines. The two equivalent engines are 
then placed at: 
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^ X GP * Y GP * 2 GP * 
avg avg avg 


^ X GP * Y GP * Z GP * 
avg avg avg 


These locations are calculated assuming equal thrust for all engines of a given 
stage, and are calculated separately for each stage. (See Figure 3-7). 


When the SRM engines are also used for moment balancing, the equivalent 
engine locations are determined by averaging the thrust-weighted gimbal point 
positions for the liquid and solid i^eled engines. The thrust weighting is 
accomplished using the first thrust value of the SRM thrust table and the liquid 
engine thrust for the first thrust event. The equivalent engine locations are 
then given by: 


* FLBS 


* SRMTHR 


and similarly for Y Gp and Z Gp 

*&vg avg 

The components of thrust at each equivalent engine location required to 
balance moments are determined using the following equations: 


T e T « — 
A 1 2 2 


where 




M,p is the total thrust moment vector 

is the total aerodynamic moment vector 
T ^ and are the equivalent engine total thrust values 


4 



V 




T is the total vehicle thrust at the current altitude. 
T 


The solution of the above equations is obtained as follows: 

Define the coordinate system xyz where the origin is located at the vehicle 
cer.ter-of-gravity and the two equivalent engines are contained in the yz plane 
(See Figure 3-8). 



Figure 3-8. CFTER OF GRAVITY/GIMBAL PLANE COORDINATE SYS I EM 


Define the p tch and yaw gitcbal angles for each engine as follows: 





i 


s 





I 


j 







T Z 

, -12 
6 - tan - — 

*2 l X 2 


The moment equations about each axis can be written: 


(T v tan 5 + T y tan 6 y )dy + (T - T y )dz = - M 

Y 1 Y 1 Y 2 Y 2 Y 1 Y 2 ** 


(T tan <5 - T tan 6 )dz - (T tan 6 + T tan 6 )dx = - M 

* X *1 *2 2 1 1 2 2 A 


(T + T )dx + (T tan 6 + T tan 6 )dy - - M 

1 Y 2 1 1 2 2 


where : 


T , T , T , T , T y , and T are the thrust components of each 
X 1 Y 1 1 X 2 2 2 

engine in the xyz coordinate system. 


^AX* ^AY* ant * **AZ are t * ie coln P onent:s the total aerodynamic moment 
in the xyz system. 


d * 5 'GPj - X CC ■ V 2 - X CG ■ 0 • 

dy : Y GPj y cc “ y gp 2 ' y cg ’ 


and 


dZ " Z GP 2 “ Z CG ^ Z GP X Z CG* * 


Since there are only three equations to determine the four gimbal angles, 
the assumption is made that 


6 Y ° 6 Y 
1 2 
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We now have four equations and four unknowns. The solution of these 
equations is presented in Appendix C and yields: 

t t r tan S tan S i 

TXX * - / ■ ' ■ - --- + , , TT 1 

/l + tan^ 6 ' + tan^ 6 ' /l + tan^ 5p t + tan^ <‘>yj 


2 


TYY = ■— 


1 

A + tar»2 6 ’ + 


— 1 1 

tan^ 6* /l + tan^ 6 1 + tan^ 6 f j 

V P_ V — ‘ 


TZZ = TYY tan 6^ 

SRH Booster Moment Balance , When the SRM booster shuttle is used, the moment 

balance equations are the same, but the values of the moments to be balanced 

and of the thrust used for moment balancing may change. With the fixed SRM 

option (M0MBAL =1), the moment produced by the SRM 1 8 is added to m „ and the 

AZ 

x f and y* thrust components produced by the SRM f s are added to the main engine 
x' and y* thrust components. With SRM’s tracking the c.g. (M0MBAL at 2), no mo- 
ment is produced by the SRM's. The x f and y f thrust components are added to 
the main engine x* and y* components (which along with the z' component, bal- 
ance the aerodynamic moments) as before. When the SRM ? s are used in conjunc- 
tion with the main engines for moment balance, the two engine equivalent 
thrust model represents the total vehicle thrust available and moments are 
balanced using this total thrust. 

Coordinated Turn Moment Balance . The thrust components Txx and Tyy are found 
from a simplified moment balance when coordinated turns are used. Since only 
pitch plane moments are encountered with this option, the moment balance pro- 
blem reduces simply to finding the angle 6 (see Figure 3-9) required to pro- 
duce no net moment on the vehicle. 

From Figure 3-9 we see that the magnitude of the moment produced by the thrust 
vector is simply 


*r“ V • 


and 
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Figure 3-9. COORDINATED TURN MOMENT BALANCE 


sin (6 - 0) • 


/ A 2 " " 2 
dx + dy 


where 


dx - x x 


d y 5 y . 

p s 


cosQ = 


/ 2 " 2 

dx + dy 


The thrust components are 


Txx = Tsin6 


Tyy = Tcos6 


and 6 is given by 


6 = 0 + (6 - 9) 
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dy 


' = cos 


— + sin 


-1 


V 1 


/ dx 2 + dy 2 


^.22 

dx + dy 


But we desire * -M^ so that 


6 = cos 


-1 


dy 


2 2 
dx + dy 


+ sin 


-1 


M 


AZ 


2 2 
T dx + dy 


the vehicle attitude angles chi-pitch, Xp* and chi-yaw, Xy* 


3.2.7 Control Low 

The control law used by RAGMOP is a set of time-dependent polynomials for 

The attitude 

polynomials are separated for each stage, so that altogether four attitude 
polynomials exist in the program - one for Xp In each stage and one for Xy 
in each stage* The current version of the program allows a quartic polynomial 
for the first stage Xp» and a quadratic polynomial for the remaining three 
control angles. The Xp profile may be discontinuous at staging, at the user's 


option, however, the Xy profile may not. 


The Xp and Xy angles, defined in 


Paragraph 3. 2.2. 5, are determined at any given time T from the polynomial 
equations : 


Xp - Xp + AP(1,1)(T-TLIFT) + AP(2,1)(T-TLIFT) 2 

BOOSTER T=TLIFT 

+ AP ( 3, 1) (T-TLIFT) 3 + AP(4,1) (T-TLIFT) 4 

where x p is the value of x p at the end of the vertical rise from the 

P T*TLIFT 
launch pad, 

and AP(1,1) through AP(4,1) are the coefficients of the polynomial, some of 
which may be zero if a lower order polynomial is desired. 


*P 

ORBITER 
where x p 

T=TSTG 


- X + AP(1,2) (T-TSTG) + AP (2, 2) (T-TSTG) 2 

P T«TSTG 

is the value of x p at the beginning of the second stage. This 

may or may not be the value of Xp at the end of the first stage, 
at the user’s option. 
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"<3t 


x y 


BOOSTER 

ORBITER 


X„ + AY(1,1)(T-TLIFT) + AY(2,i) (T-TLIFT) 2 

y T-TLIFT 

x + AY(1,2)(T-TSTG) + AY(2,2) (T-TSTG) 2 

y T=TSTG 


The polynomial coefficients used in the above equations are determined 
by curve-fitting the appropriate number of points (three for a quadratic, four 
for a cubic, etc) of \ and time. The time points are not all required to lie 
within the time spans of the particular stages; however, an even distribution 
over each stage's time span would probably be most efficient in the optimization 
scheme. The actual attitude control parameters optimized by the program 
are the Xp and Xy values at the various time points which are used to form 
the polynomials. 

3. 2. 7.1 Optional Angle-of-Attack Control for First Stage . In addition to 
the optimized x p polynomial, the vehicle pitch attitude may be controlled 
with any of three angle-of-attack options: (1) zero angle-of-attack, (2) 

zero normal force angle-of-attack, or (3) Mach number-dependent .ngle-of- 
attack. 

is found by noting that 


= 0 


°N 

o 



< i 


The Mach number dependent angle-of-attack is found using a spline 
interpolation scheme in an input Mach number/angle-of-attack table. 

Information concerning the use of these options may be found in Section 
IX, INPUT/OUTPUT. 


The zero normal force angle-of-attack 


F AFRO = 1/2p V R S(C N + C N a) 

^NORMAL R o N 


a 


AERO Hn _„. T ~N + S ° 

NORMAL o a 


and that F.^„„ = 0 when C„, + 

or, solving for a: 
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3* 2, 7. 2 Coordinated Turn Option * The coordinated turn option incorporated 
into RAGMOP takes advantage of the fact that simple relationships exist for 
describing the total aerodynamic and thrust forces on the vehicle during 
coordinated turns. The simplification occurs by noting that the relative 
velocity vector, V , lies in the vehicle x* - y' plane as shown in Figure 3-10. 

K i * 

The forces on the vehicle will also lie entirely in the x -y plane since aero- 
dynamic side forces, yawing moments, and rolling moments are all zero for zero 
sideslip angle in RAGMOP. Thus, no net aerodynamic or thrust forces on the 
vehicle are produced or required normal to the x f - y # plane. 


T + T 

yy ys 



Figure 3-10. VEHICLE PITCH PLANE 

The Xp and Xy polynomials orient the vehicle longitudinal (y f ) axis in 

the inertial reference frame. The current state determines the orientation of 

► 

the relative velocity vector, V . For coordinated turns, the plane containing 

Hs* K 

y’ and V must also contain the vertical axis of the vehicle, x'. Two possi- 
R 

bilities exist for doing this: one with positive angle-of-attack and one with 

negative angle-of-attack. Note that the pcsition of the y' axis does not change 
for the two. Rather, the vehicle is rolled about y* and V is either above 
(negative a) or below (positive a) the cockpit. We define the unit vectors 
y f and in the directions of the body y* axis and the relative velocity vec- 
tor, respectively. These are: 


y 1 = sinXpCosXy i + cosx^cosx^ j + sinXy k 


and V = i + J + v“ k * 

R R R 


■ ■ 




Consider the nonorthogona L coordinate system termed hy the unit vertois 

y* and V in the y' and V directions. The aerodynamic and thrust forces in 
K K 

this system are described as foiLows, taking note of Figure 3-11. 



The x r force (Txx + Txs - F^) may be formed by adding a torce 9) in the* 

V direction and a force in the y f direction as seen by the para I 1 eJ ogram 

K 

of Figure 3-11 where 

l l * [(Txx + Txs - F AN )/sina] V R 


l 


2 


-fc^cosa = [ (F^ - Txx - Txs) cosa/sina] y 


anc u - +cos (V • y) where the sign is specified by input. 
The total force acting on the vehicle becomes then: 


F - [Tyy t Tys - + (F^ - Txx + Txs)/tana] y + [(Txx + Txs - F AM )/sirv»]V 


AN' 


Since y and V are described in the launch plumbline XYZ system as 
R 


y “ X sinx p cosx y + Y cosx p cosx y + Z sinx y 


• . . „ W , „ U , „ V 

and V R = x v Y v Z V 
R R R 


the total force becomes 


3-J2 






t \ 


r > 


' w ^ 

F 

X 


siny cosy 
p y 


V R 





u 

F 

y 

it 

-A- 

COSY COSY 

p y 


V R 

F 


sin K ! 


V 

\ y 


y J 


K J 


The event may arise that a is too large and must be reduced in order to 

meet structural load (qa) requirements. If this is required, a rotation of the 

above y axis through the angle ( a new ~ a) is necessary, resulting in a new 

v and x • This rotation, contained in the x - y - V plane, is also shown 
P y K 

in Figure 3-12. 


By inspection of Figure 3-13, we see that 



Figure 3-12. ANGLE-OF-ATTACK LIMITING WITH COORDINATED TURNS 
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The total force acting on the vehicle in the inertial plumbline XYZ system is 
then 


F Y + F„ V„ 
y new V„ R 
■'new R 


where 

F * Tyy + Tys - + (F - Txx + T*s)/tano 

y new 


and 



(Txx + Txs - F )/sina 
AN 


The thrust and aerodynamic force components used in these equations are 
recalculated based on the new angle-of-attack. 


The algebraic sign of the angle-of-attack will be taken as a constant for 
each thrust event. A new input variable, NC00RD(15) allows the use of coor- 
dinated turns with positive or negative angle-of-attack for NC00RD (ITHR) - +1 
or -1 , respectively, where ITHR denotes the number of the thrust event. If 
NC00RD (ITHR) - 0, or if the vertical rise is not yet completed, coordinated 
turns are not used except by coincidence. 


When the coordinated turn option is used, the inertial roll attitude angle, 

/ , must be determined for print output in subroutine APRTN. Also, when the 
R 

angle-of-attack must be changed, the x p and x y angles are changed. Although 
the new values of x and X y are not required in the equations of motion, the 
values are required fir the trajectory block output of subroutine APRTN. Sub- 
routine ADER calculates the components of the body axis y system and stores them 
in YHAT(l-3) so that (see Figure 3-13): 

A A ^ A 

y = YHAT(l) i + YHAT(2) j + YHAT(3)k. 

The x and X angles are then given by 

p y 


X - tan 
P 


-1 


YHAT(l) 

YHAT(2) 
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3.3 OPTIMIZATION/RESTORATION 

RAGMOP lines a Henrrh-;ir.colorated gradient project Lon technique to com- 
pute* tin* optimum values ol a chosen net of parameters used In controlling a 
rocket ascent trajectory. As with any other s toepes t-descent method, an 
initial guess parameter set is used to start the program. Changes to this 
parameter set are calculated to meet the required end conditions and optimize 
(extremize) the payoff quantity. Since the control is limited to a parametric 
form, the solution obtained is limited to the best control of that chosen form. 
RAGMOP requires parameterization of the chi-pitch and chi-yaw attitude histories 
(see Paragraph 3.2.6) with polynomials in time, which have generally been seen 
to be in close agreement with variational solutions for the attitude histories. 

The initial guess set of parameters used to start the iteration process 
is arbitrary; however, the better the initial guess, the faster the solution 
will be found. Updating of the parameter set is accomplished in RAGMOP in two 
ways: (1) a constraint-error-reduction, or restoration, step; and (2) a payoff- 

extremizing, or optimization, step. 

3.3.1 Restorotion Step 

The first task to be accomplished by the parameter update scheme is meeting 
the desired end conations (constraints). To this end, the initial guess 
parameter set is updated using the relation 

{dB} * - H X* (X* T H X*)' 1 it, (3-1) 

where 

{ dS } is the vector of parameter update values, 

H is a diagonal weighting matrix which attempts to normalize 
the effects of the various parameters (see Paragraph 3.3.5), 

X Is the m x n matrix of partial derivatives of the constraints 
with respect to the parameters, Bj(j°l,n), 
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and 


UiT \l> 

Xg is the transpose of Xg , 


t|> is the vector of constraint errors. 

Since more parameters than constraints (n > m) are required if optimization is 
to be performed, more than one solution exists to the problem of updating the 
parameters to satisfy the constraints. The solution obtained from equation 
(3-1) is the one which minimizes the weighted sum of the squares of the param- 
eter changes, i.e., 

(d$ T H dB) is minimized. 

The minimum weighted parameter change is desirable since the update is made 
using linear prediction, and the actual problem is in general very nonlinear. 
(Small changes will be more in the region of linear prediction) . 


The restoration step of equation (3-1) is scaled whenever any of the 
constraint errors is greater in absolute value than fifty times the acceptable 
error for a converged run: i.e. 

dB - E 2 (-HXg(Xg T HXg) -1 t) , (3-2) 


where E 2 is the restoration step scale factor. The value of is found 
assuming a quadratic relationship between a performance index P^ and the scale 
factor E^: 


+ aE„ + bE„ 




where 


P 


is a constraint performance index for the total constraint errors. 



and 


a and b are the coefficients in the quadratic relationship between 

P and E 0 . 

2 


The constraint performance index used in RAGMOP is 


m 


l k ♦ 


i»l 


i i 


where 

|gnu ± I 

k i " END^ 

{GNU} = - (A* T HA*) (Af HAjjjr 1 , 

END^ is the maximum acceptable absolute arror in the i ^ constraint 
for a converged run, 

<j) 

is the vector of partial derivatives of the payoff with respect 
to the parameters, 

and 

<t> is the payoff. 

The assumed quadratic relationship between P^ and E^ is established using 
three bits of information: (1) the value of P^ when 35 0, (2) the slope of 

the P^(E 2 ) function at E 2 ® 0, and (3) the value of at some value of 

E 2 ^ 0. The first two bits of information are readily available since 

P ■ P = P 

tfr ^ \h \b 

E,“0 ^current v o 


and 
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k-° 


since 


/ GNU . > 

2 U END? V 


i / e 2 -o 


(*l\ 

I 3E J 
\ 'E o «0 


The third bit of information is the value of P, at some nonzero value of E 0 . 

4 » 2 

This is obtained by choosing E«, = .1, using the parameter update with that 

value of E 0 , running the trajectory, and evaluating P. . The quadratic rela- 

1 *E 

tion between F. and E_ is, then: 2 

* 2 

3P * 2 
P, * p , + rp- E, + bE_ 

i'c «> 3E 9 2 2 

Ej o 2 


which yields 


P. ~ P. - 2P. K„ P ( - P, (1 + 2K ) 


E 2 


,: 2 ^ 


Since we desire all of the constraint errors to be nulled we have P, 

i 

so that: 


0 * P, * 2P. E 0 + bE 0 
^ ^22 
c o 


desired 


or solving for E 
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E 


2 


E 


2 



A 


4P. 


2b 


(P 


2 


4bP, 


b) 


b 


The value of E 2 used in RAGMOP is based on the value of b. For b < 0, 

E~ = P, + /P (P, - b) 

2 

o o o 

b 

and for b > 0 E« = ^o 

b 

The value of E 2 is limited by RAGMOP to .05 < E 2 < 1.0. 

3.3.2 Optimization Step 

The second task to be accomplished by the parameter update scheme (after 
the constraints are satisfied) is the ext realization of some payoff quantity - 
The step which accomplishes this is given by 

dS ■ QY {-„(*« - X*(** T HA*)' 1 »« HA*))} (3-3) 


where 


QY is a step scale factor to be determined by a search. 

This step assumes linearity of the payoff and constraints with respect to the 
parameters and also assumes that the optimum on the intersection of the con- 
straint hyperplanes (which describes the allowable surface of parameter values 
over which the optimum payoff is sought) can be approached using a linear step 
Since the linearity assumption is false, the step is scaled to Increase the 
payoff as much as pospible without excessive increase in the constraint errors 
In other words, the attempt is made to remain within some defined region of 
the constraint boundary during the optimization step by limiting the step size 
taken • 
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The value of the optimization step scale factor QY is found using a one- 
dimensional search on the change in value of a composite payoff index, A$: 

A$ - A$ + 


where 


A$ 3 $ — $ 


QY=0 


AiJ/ = ^ ^ 


QY s O 


and 


v * -(Xg T (Xg T H *g) 1 GNU of Paragraph 3.3.1). 

The search begins with a value of QY calculated to give some minimum parameter 
change, below which the significance of the search is lost due to machine noise 
and/or roundoff error. (A change in liftoff weight of .1 lb for a 3 million- 
pound vehicle is not significant, for example.) The parameter update is tem- 
porarily performed with this QY, a trajectory run is made, and the composite 
payoff index evaluated. The value of QY is then increased, the temporary 
parameter update performed, the trajectory run made, and the composite payoff 
index reevaluated. This procedure continues until the value of the payoff 
index increases, denoting that the minimum has been passed. If more than four 
values of QY are tried before the minimum is passed, only the most recent 
four points are used in a cubic curve-fit to locate the minimum, if the 
minimum is passed before four points are obtained, three points are run and a 
quadratic curve-fit performed for the minimum. The result (QY) from this curve- 
fit is then used to form a fourth data point for a cubic curve fit. Once the 
first cubic curve-fit is performed, the resultant QY is used for another data 
point, and a second cubic curve-fit is performed for the value of QY that 
minimizes the payoff index $ (Note that for maximized payoff a negative sign 
is introduced so that the minimization logic is still valid). 
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3. 3. 2.1 Cubic Curve Fit 

The cubic form is written 


Y = a X 3 + bX 2 + cX 4 d, where Y = $ and X = QY. 


(3-4) 


The value d is just the value of Y at some base point, (usually at X = 0, 
but not necessarily so) so choosing the fourth data point as the base point 
we have: 


Y - Y 4 = a(X - X 4 ) 3 + b(X - X 4 ) 2 + c(X - X 4 > 


(3-5) 


At the data points (X^ Y^) (X 2> Y 2 ) (X^, Y^) (X 4> Y^) we can write: 

3 2 


H - - 4 - a(X 1 - X 4 ) + bCXj^ - X 4 ) + c(X L - X 4 ) ^ 


Y, - Y. 


Y 2 - Y 4 = a(X 2 - X 4 ) 3 + b(X 2 - X 4 ) 2 + c(X 2 - X,) 
y 3 - Y 4 = a(X 3 - X 4 ) 3 + b(X 3 - X 4 ) 2 + c(X 3 - X 4 ). . 


(3-6) 


Ue now have three equations in the three unknown coefficients a,b, and c. 
The solution for these coefficients follows: 


dXi 

III 

X 

1 

X 

dY 

dX 2 

5 X 2 - X 4 

dY 

dx 3 

. x 3 - x 4 

dY 


1 *1 4 

s y - Y 

2 *2 4 

= y - Y 

3 3 4 


then Equation (3-6) can be written: 


dY 


dX 


1 = adX 2 + bdX, + c 


(3-7) 


dX. 


adX 2 + bdX 2 + c 
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fa 

dX. 


= adX 3 + bdX 3 + c 


(3-9) 


Subtracting Equation (3-8) from Equation (3-7) gives 


dY dY 

dx^ * d5£ = a(dx i ' " 2 > + b(dx i - dx 2> 


(3-10) 


anc similarly (3-9) from (3-7) gives 


dY l dY 3 2 2 

dx7 - dx! = a(dX l - dX 3 > + b (dX l - 


*1 3 

which reduce to 


(3-11) 


dY, 

dY. 

1 

2 

dX, 

dX_ 

1 

2 

dY, 

dY 0 

1 

3 


dX x dX 3 


) / (dx 3 - dx 2 ) 


) / (dx x - dx 3 ) 


a(dX^ + dX^) + b 


a(dX A + dX 3 ) + b 


(3-12) 


(3-13) 


Defining dX^ 5 dX 3 - dX 2 


and dX 13 = dX 3 - dXj 


and subtracting (3-13) from (3-12) we have 


d* dY 

( ) 

1 dX 3 dX 2 ’ 


a(dX 2 - dX 3 ) 


dX 


12 


fl fa 

( dX 3 “ dX ;} ) 


dX 


13 


(3-14) 


and now defining dX 23 " *^2 ” dX 3 
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we have 


dY, 


dx i 

a = 1 


dY <IY 

dX') 7 dX J2 -<5^ 


dY 

dxl> 7 dx i3 


dX 


23 


( 3 - 15 ) 


and then from (3-13) 


dY l dY 3 

b = ( dX- - s;> 7 "13 - a(dX l + dx 3 ) 


(3-16) 


and from (3-9) 
dY 3 

c = + adX 2 + bdX 

dX 3 3 J 


(3-17) 


Using (3-15), (3-16) and (3-17) in (3-5) we have 


or 


^ILN 


^SlIN 


n. 

-2b + 4b - 12ac 


6a 

-b -K b 2 - 3ac 


3a 


(3-18) 


Thus (16) yields the extrema of the function. In order to choose L 'tween 
the + or - sign we look at the second derivative at the extrema: 



6aX + 2b = 


/— 2 

-6ab 4- 6a b — 3ac 


3a 


which becomes 


+ 2b > 0 
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from which 


> 0 

- 3ac 

the choice of 


QY 


min 


3a 


the + sign is obvious and 


3a c 


This value of QY is then used to perform a permanent parameter update (the 
optimization step). The trajectory is run with the permanent update, a new 
set of partial derivatives is obtained at the new control state, and the 
convergence test is made. If the solution has not yet converged, another QY 
search will be performed as soon as the constraint errors from the optimization 
step are restored to acceptable levels (if they ere not already acceptable). 

The search is terminated prior to pessing the minimum if either of two 
criteria are not satisfied: (1) parameter changes must be below some maximum 

level, and (2) the constraint errors nust not be excessive. The parameter 
change is monitored by determining a maximum acceptable change for a single 
parameter which directly affects the payoff. The constraint errors are moni- 
tored by determining the maximum act 'ptable error of a single constraint. If 
the parameter change or the constraint e. ror exceeds the maximum value, the 
search is terminated. The current value of QY is then used for the optimization 
step calculation rather than a curve-fitted value. 

3.3.4 Porti a I Derivatives 

The partial derivatives used in the parameter update equations 


•P 

r 



n 2 

3^ 

m 


3e x * ’ ' 

~ae 1 

a »i 






3^ 



c 

QQ 
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fii- 1 


38, 


38„ 


3(fr 


36 


n 


are determined using numerical difference techniques. Each parameter is 
perturbed by some slight amount, a trajectory is run, and the effect of the 
perturbation on the equality constraints and on the payoff is calculated. This 
process must be performed separately for each parameter. The partial deriv- 
atives are then found as the ratio of the change in the constraint or payoff 
to the change in the parameter: 


alT 




- 4>4 


8.+A8. 
-J L 


8, 


A8. 


or 



V+A8. V 
1 . .1 i 

A6. 

J 


Derivatives formed in this manner (forward differences) are used whenever a 
restoration step is made. If the optimization step (QY) search is to be 
performed (signalled by all constraint errors being within the acceptable limits 
for a converged run), the parameters are also perturbed in the negative direc- 
tion and the average slope (partial derivative) is found: 


if) - ijj 

i0 1 +A6 i i B 1 -AB 
W 5 ’ 2A6 j 
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and 


* V AB 1 V A8 i 


38 


j 


2AB , 


This method is used since the optimization step (QY) search is more sensitive 
to errors in the partial derivatives than the restoration step. (A consider- 
able savings in computer time results from using only the forward differences 
whenever practical.) 

The perturbation values used in RAGMOP are: 

A burn times » 1 second 

A liftoff weight * .01% 

A launch azimuth = .25° 

A first stage Xp ® .1° ^ 

A first stage Xy ® *5° 

A second stage Xp * *1° 

A second stage Xy ■ *5°^ 

The one exception to the above method of determining the partial deriva- 
tives occurs when the final thrust event bum time is optimized without the 
payload option. In this case, the partial derivatives are obtained analytically 
rather chan by a trajectory integration. The approach used is that 


Preset to these values. May also 
be changed by input. 


KM) 

3 (TAUT) 


KMl 

3 x 


3 x 


3 (TAUT) 


where 

x denotes the final state (position, velocity mass). 
The derivatives 

3 x 
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4 


are presented in Appendix G. The use of this formulation both saves computer 
time and increases the accuracy of the partial derivatives by eliminating the 
numerical difference technique. 

3.3.5 The Weighting Matrix 

The purpose of the weighting matrix is to normalize the partial derivatives 
of the constraints with respect to the various parameters. The numerical 
value of each partial derivative 


^ . 



is dependent upon the units chosen for the parameter 8^ . Consider for example 
the case of the partial derivative of the radius vector at orbital injection 
with respect to one of the values used to form the first stage x polynomial. 
At some time suppose the value of this partial derivative is 



1000 meters /degree 


If the units of Xp were radians rather than degrees, the value would be 17.45 
instead of 1000, and the square of the partial (used in the parameter update 
scheme) would be 304.57 instead of 10^. Now suppose that the partial derivative 
of the radius vector with respect to liftoff weight is 

t~~ r = 1 meter/kilogram 


The square of this is 1. With Xp in radians, the ratio of the square of the 
liftoff weight partial to the Xp partial is 1/304.57 = .00328, but with Xp in 
degrees this ratio is 1/10^ =» .000001. Thus, merely by changing the units of 
Xp, the effect of the x partial derivative on the sum of the squares 

2 2 
j)R_ 3R 

3X 9W0l 

P 

is changed from the third significant digit to the sixth. It is easily seen 
that for some parameters and constraints, the significance of one or more 
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parameters on a particular constraint could be completely lost on an eight 
digit computer when summing the squares as shown above. The following para- 
graph explains the significance of this sum. 


The matrix 




HA 


♦ 

8 


must be inverted in the parameter update calculation for both the restoration 
step and the optimization step, /lie diagonal elements of this matrix contain 
the weighted sum of the squares of the partial derivatives of the parameters 
with respect to each constraint, i.e.. 


3i|> 


3 if 


“ H 11 + H 22 38 


+ H 


nn 


3^ 

38 


Without a weighting matrix the numerical value of the terms comprising each 




are dependent upon two factors: (1) the actual effect of the various 


MJ th 

parameters on the j constraint, and (2) the units chosen for the parameters. 
The weighting matrix attempts to eliminate the second factor by normalizing 
all of the squared terms to the same order of magnitude. 

The weighting matrix is calculated as follows: 

Define the vector 


38 ! 


2 \ 


V = 


3*j 


36, 


3-49 





* 


<1, 


s 




fz 



where it is assumed that all the V > 0. Using V, a second vector S is 
chosen such that 



ioe ^ 


a 1 * 1 

< 8T> 

j 


2 


V 


i 


m 


In other words, is the average base ten logarithm of the ratii of the squares 
of the constraint partial derivatives of each parameter with respect to the 
first parameter. Thus, ■ 1, 

The first parameter is arbitrarily chosen as the basis for comparison, how- 
ever, any parameter with all the > 0 could be used. (The parameter chosen 
affects the scaling of QY in the optimization step.) 


The weighting matrix U is then given by 


"jj ' V 10 3 


where is an additional weighting factor found by experience to improve the 
performance of certain parameters.. This definition of the weighting matrix 
has been seen to contribute to essentially uniform convergence of the para- 
meters. 


The weighting matrix calculation is performed at each restoration step 
after the initial guess (nominal) trajectory until an optimization step scale 
factor (QY) search is completed. At this time the weighting matrix is updated 
using the Fletcher and Powell form of Davidon*s update formula such that: 
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where 


A*, 


- D4> 

" D6 i 


D<t> 


DB 


i-1 


and 


AB = Bi - B ± _ 1 


^ * total derivative of the payoff with respect to parameters. 

Using the updated weighting matrix, the restoration step is again per- 
formed unt?l the constraints are satisfied. At this time, the convergence 
tests as described in subsection 3.3.6 are tested, if the tests ar met con- 
vergence is denoted and a final trajectory printed. If the tests are u t met 
another optimization step search is begun and the weighting matrix updated 
upon successful completion of the search. This process of updating the 
weighting matrix is continued until either the convergence tests ar -net or 
th i updating process is performed more times than the difference between the 
number of parameters and the number constraints. In the latter case the 
weighting matrix is reinitialized based upon the aforementioned procedure. 
Experience has shown that the reinitialization process speeds the convergence 
since this process tends to eliminate the inherent numerical errors. 

3.3.6 Convergence Test 

Three criteria must be met at the end of a QY search (See Paragraph 3.3.2) 
in order for a given trajectory solution to be considered converged: (1) the 

constraint errors must all be within acceptable limits, (2) the scaled total 
derivative of the payoff with respect to each parameter ( p con ) must be 

acceptably small, (3) the parameter changes computed using the curve-fit value 
of QY must be trivial. 

The acceptable limits on the constraint errors are input values. 
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The total derivatives of the payoff with respect to each of the parameter 
are given by: 
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The convergence parameters mentioned in (2) above are given by 


P 

coi^ 


D<j) 



P scale i 


where P sca l e 1 ls ‘ 


calculating 


M. 


absolute value of largest of the terms added together in 
This scaling assures that each is zero to the same 


number of significant digits. The test on P in RAGMOP is set at 

Dd> Con i 

P „ < .025. (On the optimum = 0 and P ■ 0.) 

con r D8. con. 


The absolute values of the parameter changes at the end of a QY search 
(i.e. using the value of QY produced by the search) must all be trivial, 
where trivial for this purpose is defined as: 


A burn times < IX 

A liftoff weight < .001% 

A launch azimuth < .00085' 


A y values 
P 

A v values 

y 


< IX 

< 0 . 1 % 


3.4 INTEGRATION PACKAGE 

Numerical integration is performed in RAGMOP using the DESOLV integration 
package The DESOLV package provides executive control of the integration 
of the equations of motion from initialization to termination. Three integra- 
tion schemes are available in the DESOLV package: 1) a fixed step fourth 
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order Runge-Kutta, 2) a variable step-size, variable order. Adaras -Moulton 
(maximum order) with fourth order Runge-Kutta starter, ami 1) a fixed step 
sire, variable order Adams predictor (maximum H*"' 1 order) with fourth order 
Kunge Kutta starter. 

Conicux of the trajectory integration by DESOLV is performed using a 
series of integration interrupt triggers. These triggers serve two purposes: 

(1) to allow a print of the state at various times during the trajectory, and 

(2) to allow changes in the equations of motion to be made during the trajectory 
integration (such as termination of the thrust events based cn time, fuel, 
relative velocity, or acceleration) . The integration package recognizes the 
difference between these two types of triggers, and will restart the integration 
for the second type. 
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MANUAL 





Section V 

MAIN PROGRAM FLOW AND OPERATION 


This section describes the flow and operation of the RAGMOP main program. 
Individual subroutines, treated as "black box" units in this section, are 
described in detail in Section VI. The main program is best described in six 
segments : 

• Initialization, 

• Trajectory integration and evaluation, 

• Convergence test , 

• Partial derivative calculation, 

• Parameter update, 

• Plot and table output. 

These six segments are described In detail in the following paragraphs. A 
flow diagram of the main program is presented in Figure 5-1. 

5.1 INITIALIZATION 

The main program begins by presetting a number of variables and constants 
used in the program. Section IX presents the preset values of all input 
variables. In addition to the preset input variables, several program control 
variables are initialized and the weighting matrix is preset to an n x n 
identity matrix where n is the number of parameters to be optimized. The input 
subroutine AINIT is then called to read the input data for the case to be run. 
This input data Includes a description of the vehicle, the type of trajectory 
to be flown, the parameters to be optimized, and the output desired (see Sec- 
tion IX). The maximum allowable parameter change during the optimization step 
is calculated based on which parameters are used. 

5.2 TRAJECTORY INTEGRATION AND EVALUATION 

The trajectory integration and evaluation is performed by calling AF0RUN 
and AF0RND. Upon completion of these two subroutines, the values of the 
c*..istraint ertors and the payoff are known for a particular parameter se.. The 
constraint errors are tested to determine whether any exceed the acceptable 
level in END and fifty times that level. 
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Figure 5-1. RAGMOP MAIN PROGRAM FLOW DIAGRAM 
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PARAMETER UPDATE (Continue) 



Figure 5-1. RAGMOP MAIN PROGRAM FLOW DIAGRAM (Concluded) 
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5.3 CONVERGENCE TEST 

A test for convergence is the next step. A run is considered converged 
when, using the parameter update at the completion of an optimization step 
scale-factor search, the constraint errors are small and the convergence tests 
of subroutine ANEVCH have passed (see Sections III and VI). If the convergence 
test is passed, NMAX is set to zero and the trajectory integration is performed 
again to obtain a detailed print of the final trajectory and to create a 
trajectory tape for the plot and table routine if they are to be used. If the 
convergence test is not passed, the payoff and the constraint error values are 
stored in PHITES for later use as baseline values in the restoration and opti- 
mization step scaling. 

5.4 PARTIAL DERIVATIVE CALCULATION 

The ^rtial derivatives of the constraint errors and the payoff with 
respect to the parameters are then determined by calling subroutine BADLX. 

Upon the return to MAIN, the partial derivatives will be stored in the array 
XLAMB. These derivatives are determined using forward numerical differences 
when the constraint errors are large and central differences when the con- 
straints ar* within acceptable convergence limits 'e.g. at the beginning of 
an optimization step scale-factor search). 

5.5 PARAMETER UPDATE 

The parameter update is performed in one of two ways: (1) as a constraint- 
restoring step, or (2) as an optimization step after the constraints are met. 


Once the partials have been obtained, the values of the parameters used 
for the last trajectory (before the partial derivative trajectories in BADLX) 
are stored in the vector PARSAV. 

The variable NTNTST is checked to determint whether any of the constraints 
were greater in absolute value than fifty times the allowable limits in END 
(indicated by NTNTST»1). If NTNTST*!, the restoration step is scaled to 
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minimize the total constraint error based on a quadratic assumption for the 
relationship between a constraint performance Index and the scale factor DP2. 

If NTNTST^l (indicating that ail constraint errors are less than 50 times the 
acceptable limits) the restoration step is not scaled. For either case sub- 
routine ANEWCH is called to calculate the matrix products and inverses required 
in the parameter update. When the restoration step is not scaled and the 
constraint errors are n. yet sufficiently small to perform an optimization 
step search , this call to ANEWCH (with DP2®1.) is also used to perform a 
permanent parameter update and the flow returns to the trajectory integration 
of Paragraph 5.2. 

When scaling of the restoration step is required , subroutine ANEWCH is 
called with DP2=.l, a temporary parameter update is performed, and the trajec- 
tory is run in order to obtain the value of a constraint performance index. 

This value serves as a data point in calculating the value of the scale factor 
(DP2) that minimizes the performance index for use in the permanent parameter 
update. A quadratic relationship is assumed between the performance index and 
the scale factor for the purpose of this scaling. After performing the 
parameter update with the calculated scale factor, the flow returns to the 
trajectory integration of Paragraph 5.2. 

If all constraint errors are less than the acceptable levels for a con- 
verged run, the optimization step scale-factor (QY) search is begun. The value 
of QY which will produce a minimum acceptable parameter change is found and 
used as the starting value for the search. This is done to eliminate searching 
in the machine "noise" region with very small parameter changes which might 
result from an arbitrarily small starting value for QY. The parameter update 
is performed and a composite payoff index is checked for improvement over the 
baseline value at the beginning of the search. When the minimum of the payoff 
index is passed, and at least three values of QY have been tried, a curve fit 
scheme predicts the value of QY which yields the minimum. If the minimum has 
been passed with the third data point, a quadratic fit for QY is made. If 
after three values of QY have been tried, the payoff index is still decreasing, 
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the value of QY is repeatedly increased until an unacceptably large parameter 
change or constraint error is obtained or until the minimum is passed, with 
the most recent four data points (QY and payoff index) retained . _ each step. 

When the payoff index increases, (indicating that the minimum has been 
passed) or after the quadratic fit mentioned above when the minimum three 
data points are used to predict QY (which is then used to produce a fourth 
point), a cubic polynomial is used to estimate the peak value of QY. This 
value is then used to produce another data point, the first point is dropped, 
and the cubic fit repeated one time. The parameter update scheme is then used 
to permanently change the parameter set with the final value of QY and flow 
returns to the trajectory integration of Paragraph 5.2. The value of QY in 
subsequent restoration steps will be zero. The search is terminated early 
(before the minimum is found) if the test parameter change (W01 or the last 
optimized TAUT) becomes too large or if the error in the radius vector at 
insertion is greater *han 1000 meters. 

5.6 PLOT AND TABLE OUTPUT 

Once a converged run is obtained, or when the maximum number of iterations 
has been reached, a test is made of the input variable NTABLE to determine 
whether special output tables (NTABLE-1) or plots (NTABLE=2) or both (NTABLE=3) 
are desired. If tables are desired, subroutine BOPTBL is called. If plots 
are desired, subroutine BOPPLT is called. The Input variable LAST is then 
decremented (LAST^LAST-l) and checked for multiple cases (LAST>0). If another 
case is to be run, the flow returns to the initialization of Paragraph 5.1; 
otherwise the program stops. 
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Section VI 

SUBROUTINE DESCRIPTIONS 


This section describes in detail the subroutines (with exception of the 
integration package) used in the RAGMOP program. For each subroutine, the 
functional flow is given as a word description and a macro-flow diagram. The 
integration package subroutines DPIR, DES0LV, and RTMRK are not presented in 
a detailed fashion since documentation for these routines is available in 
reference 17. 


6.1 ALPHABETICALLY ORDERED LIST OF SUBROUTINE NAMES AND FUNCTIONS 


Subroutine Name Function 


ACST0P Evaluates the payoff and determines the value of termi- 

nal and intermediate equality constraints. 

ADER Evaluates the state derivatives. 

ADERl Evaluates time dependent variables used in ADER. 

AF0RND Controls the call to ACST0P to determine terminal con- 

ditions, computes errors from desired values, and 
prints terminal summary. 

AF0RUN Controls the setup logic for the trajectory integration 

and calls the integration package. 

AGE0 Evaluates oblate gravitational properties. 


AINIT 

AKALT 

AMACH 

AMULG 

ANEWCH 

APRTN 

AQMAX 

ASIMP 

ATHREV 


Reads input data, performs initialization. 

Trigger routine for 10 km and 14 km altitude printout. 
Trigger routine for MACH *■ 1 printout. 

Linear interpolation scheme. 

Evaluates the optimization and restoration step equations 
to find the required changes in the controlling param- 
eters. NNEWCH is an entry point in this subroutine. 

Performs output editing, prints data, and writes output 
tape for output tables and plots. 

Trigger routine for maximum dynamic pressure printout. 

Determines the flyback range using empirical data 
interpolation. 

Thrust event control routine; sets up vehicle geometry, 
aerodynamics, number of engines, thrust limits, flow- 
rates, and trigger flags for thrust event cutoff. 
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Subroutine Name 


Function 


ATILT 


AXPRT 

BADLX 

B0PPLT 

B0PTBL 

CHIP0L 

DES0LV 

DPIR 

FIND 

GLIKT 

MATINV 

NNEWCH 

PRA63 

RTMRK 
SPLINE 
SPLINE 2 
TEST 
TRIVLP 


Trigger routine for end of vertical rise (begin tilt 
over or attitude control). The subroutines AKALT, 

AMACH, AQMAX , AXPRT, AND GLIMT, are contained within this 
subroutine as entry points. 

Trigger routine for normal printout. 

Determines the Influence coefficients (partial deriv- 
atives) for computation of the changes in the 
controlling parameters. 

Generates output plots on CALCOMP plotter. 

Generates output table summary. 

Evaluates the pitch and yaw attitude polynomials. 

♦Integration package 

Determines attitude polynomials based on attitude time 
history information. 

Trigger routine for acceleration limit printout. 

Performs matrix inversion. 

Determines parameter update step during one dimensional 
search. 

Evaluates the atmospheric parameters as well as their 
partial derivatives w.r.t. altitude using the spline 
interpolation method of subroutine SPLINE. 

Part of integration package which flags termination of 
the trajectory integration. 

Inteipolation routine using SPLINE method for aerodynamic 
coefficients. ** 

Interpolation routine using SPLINE method for SRM 
thrust and weight drop tables**. 

Linear scheme used in ASIMP for determining slope of the 
dependent variable w.r.t. the independent variable. 

Trivariant lookup (of flyback range) for ASIMP. 


6.2 ALPHABETICALLY ORDERED LIST OF SYSTEM SUBROUTINE NAMES AND FUNCTIONS 

Subroutine Name • Function 

AC0S Determines arccosine (rad). 

AINT Truncates real variables to integer form. 

AL0G Determine the logarithm of a variable (base e) . 

* Documentation of the integration package can be found in reference 17. 

**See Appendix C for further information concerning interpolation techniques . 
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Subroutine Name 


Function 


AL0G1O 


Determine the logarithm of a variable (base 10) . 


AMAX1 

Selects maximum value from argument list. 

ASIN 

Determines arcsine (rad) 

ATAN 

Determine arctangent based on one argument 
- f 1 0 < j (rad) . 

ATAN2 

Determine arctangent based on two arguments 
- ir £ 0 <_ it (rad) . 

AXIS 

CALCOMP system routine, draws axis with "tic" marks and 
scales at every inch. 

CALEND 

CALCOMP plotter routine: plot an ending identification 

frame and end the current plot file. 

CALID 

CALCOMP plotter routine: initializes the CALCOM' J plotti 

uses interface package. 

cos 

Determines the cosine of the argument. 

DC0S 

Double precision C0S. 

DSIN 

Double precision SIN. 

DSQRT 

Double precision 3QRT. 

EXP 

Exponential function. 

LINF 

CALCOMP plotter routine: plots a line or curve from a 

set of points. 

NAMELIST 

Variable format input routine that reads data cards and 
places the values in the proper common blocks. 

PL0T 

CALCOMP plotter routine: moves pen to a new coordinate 

position, or establishes a new reference position. 

SCALE 

CALCOMP plotter routine: obtains maximum and minimum 

scale values to adjust plot values. 

SIN 

Determines the sine of the argument. 

SQRT 

Determines the square root of the argument. 

SYMBOL 

CALCOMP plotter routine: draws a symbolic label at a 

specified position. 




is 


6.3 DETAILED DESCRIPTION OP SUBROUTINE 

The following paragraphs describe each subroutine in detail. 


6-3 


f 



6.3.1 ACST0P 

Subroutine Identification 

• Title 
ACST0P 

• Calling sequence 

CALL ACST0P (M, VALE, I), where M - option flag see Table 6-1 

VALE - value of option used 
I - not used. 

Function 

Evaluates the payoff and determines the values of the terminal and 
intermediate equality constraints. 

ACST0P also provides the analytical partial derivatives of the constraints, 
with respect to the state variables when fin-’ burn time is a parameter and 
payload (as opposed to final mass) Is Hit . > ■ . i it or payoff. 

Functional Flow 

Called by the subroutine AF0RND at the -nd f ■ trajectory integration 
to evaluate the current payoff and to do ■' cmrat ra'i of the tt'.inai and 
intermediate constraints. The calling argument, M, of the subroutine denotes 
the calculation require^ of the routine, and V. LE denotes the value produced 
by that calculation. The options for M are defined in Tabl •>!. 


Table 6-1. FUNCTION LIBRARY FOR CONSTRAINTS 


Code No. (M) 

Function 

1 

Payload (or final mass) 

2 

Inertial velocity 

3 

Inertial flight path angle 

4 

Radius 

5 

Energy 

6 

Angular Momentum 

7 

Inertial longitude 

8 

Inertial heading angle 

9 

Colatitude 

10 

Inclination 

11 

Line of nodes 

12 

Semi 1 ; us rectum 

13 

Eccentricity 
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Table 6-1. (Concluded) 


14 

Bum Time 

15 

Maximum dynamic pressure 

16 

True anomaly 

17 

Argument of perigee 

18 

Not used 

19 

Not used 

20 

Flyback range 


The value of M is determined by input of the integer variables KCDPHI 
and KCDRES. 


ACtTOP also proviues the analytical partial derivatives of the constraints 
with respect to the state variables when final burn time is a parametric and 
pcvload (as opposed to final mass) is not a constraint or payoff. 


Functional Flow Diagram of ACSTOP 
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4.3 2 ADER 

Subroutine Identification 

• Title 
ADER 

• Calling Sequence 
l JL ADER 


t 


Function 

Evaluates the state derivatives 
Functional Flow 

Called by the subroutines AF0RUN, ATHREV, and the integration routine 
(DP1R) to calculate the values and derivatives of variables which are not 
totally dependent cn time. These incxude the current velocity and acceleration 
of the vehicle * the mass of the vehicle, the thrust, propellant flowrate, the 
aerodynamic forces and moments, the thrust vectoring required to balance 
moments, angle-of-attack, sides li? angle, dynamic pressure -ad its deriv- 
atives, and the attitude requirements when an angle-of-attack history is 
specified or when intermediate constraints of qa or qB are exceeded. 

The first part of the subroutine finds the magnitudes of the vehicle 
position and velocity vectors. The gravitational acceleration is found by 
calling subroutine AG\0. The altitude is found using either an oblate or a 
sph ~ j cal earth model, as desired. The inertial velocity of the atmosphere 
at the current vehicle altitude due to the rotation of the earth is calculated, 
and wind bias values (if any) are added. The difference betweer the veh*r’..*s 
Inertial velocity and thaL of the atmosphere is then used to determine the 
lelallvr velocity. The atmospheric routine, PRA61, is called to determine the 
presume, density, spo«'d ot sound, and partial derivative of density w.r.t. 
altitude at the current altitude. Mach number is calculated and used for the 
interpolation of the aerodynamic force and moment coefficients. Axial force 
uue to ^ase pressvre is then determined by interpolation based on current 
altitude. Dynamic pressure, q, is determined based on the density and the 
relative velocity. Angle-of-attack and sideslip angles are determined by either 
coordinated turn requirements or transforming the components of the relative 
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velocity into the body axis system. The value of the input variable FACT is 
examined to determine whether an angle-of-attack profile is required. The 
values of the structural load indicators, qa and q8» are checked to determine 
violation of input constraints. If the variables violate the constraints, a 
and/or 8 are set e ual to the constrained values and the values of pitch and/ 
or yaw attitude are calculated to satisfy the constraint relationships, or the 
desired profile, overriding predetermined values. 

The next portion of the subroutine calculates the aerodynamic and thrust 
forces and moments acting on the vehicle. Aerodynamic force and moment coef- 
ficients are found using the zero angle-of-attack coefficients, the slope of 
coefficients w.r.t. angle-of-attack or sideslip, and the aerodynamic moment 
arms. The aerodynamic forces and moments are calculated with respect to the 
body axes. 

The value of the thrust at the current altitude is found and tested for 
order of magnitude. If the thrust is too low (arbitrarily chosen 10,000 
newtons) or no moment balance is required, the moment balance scheme is by- 
passed. If sufficient thrust exists, and moment balance i reouired, the 
thrust components needed to balance the moments are found depending on the 
type of moment balance as given by the input variable M0 HBAl.. A two- 
engine equivalent is used in the moment balance scheme for the controllable 
engines regardless of the actual number of engines. If acceleration (g) 
limiting j: required, special calculations are made to determine the thrust 
requirements and propellant flowrate of the main (orbiter) engines. Only 
liquid fuel engines are throttled, and these are throttled only to zero thrust. 

If the SRM thrust alone causes the g-limit to be exceeded no further action 
is taken. 

The total forces acting on the vehicle are calculated and transformed 
to the iner plumbline system. Total acceleration in the plumbline system 
is found by summation of the aerodynamic, thrust, and gravitational accelerations. 
The final portion of the subroutine calculates the derivatives of the velocity 
losses and of the dynamic pressure. (See Appendix F for an explanation of 
the velocity loss equations). 



t 
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6.3.3 ADERl 


Subroutine Description 

• Title 
ADER1 

• Calling Sequence 
CALL ADER1 


Function 

Evaluates time dependent variables used in ADER. 

Functional Flow 

Calculates all time dependent variables, i.e. thrust, weight, and attitude 
control values. The subroutine is called initially by AF0RUN and subsequently 
by the integration package. 'Hiis subroutine must be called at the beginning 
of each integration time step and called before ADER, since calculations 
in ADER require values computed in ADER1. 

The first step of this subroutine is to test to see if SRM engines are 
being used. If so, the value of SRM thrust and weight loss are determined 
from SPLINE interpolation of time dependent tables. The value of the pro- 
pellant flowrate for the liquid engines is stored into the variable DVAR(7) 
for integration of mass loss. A check is made to determine whether the 
liftoff phase has been terminated (KS1>2) and if so, the appropriate attitude 
values are found. Finally the location of the center-of-gravity at the 
present vehicle weight is found, cs are the distances of the engine gimbal 
points from the present center-of-gravity position and several relations 
among these distances used In the moment balance equations. 
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Functional Flow Diagram of ADER1 
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6.3.4 AF0RND 

Subroutine Identification 

• Title 
AF0RND 

• Calling Sequence 
CALL AF0RND 


4 


.1 

1 


N 


t ; 




i 




c 


t 


f 



Function 

Controls the call to ACST0P to determine terminal conditions, computes 
errors from desired values, and prints terminal s\ unary. 

Functional Flow 

AF0RND is called at the end of each trajectory run. The subroutine 
< ecks the constraints by calling the subroutine ACST0P for the payoff and 
terminal constraints requested by input and evaluates the differences between 
the desired terminal values and the actual values. If AF0RND has been called 
□y BADLX (indicated by LSB»1) the constraint errors are evaluated and returned 
to BADLX to calculate influence coefficients. If AF0RND has been called by 
MAIN (LSB=0) an additional trajectory summary printout is printed and returned 
to MAIN. 


i 



i 
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Functional Flow Diagram of AF^RND 





6.3.5 AF0RUN 

Subroutine Identification 


• Title 
AF0RUN 

• Calling Sequence 
CALL AF0RUN 


Function 

AF0RUN controls the setup logic for the trajectory integration. 

Functional Flow 

Calling the AF0RUN subroutine is synonymous to calling the integration 
of the trajectory from the initial conditions to termination. The subroutine 
is called either from MAIN or from BADLX. 

The AF0RUN subroutine first calculates a transformation matrix for 
rotation from an inertial equatorial coordinate system to the inertial plumb- 
line system. If the trajectory is initiated at liftoff (JUMP=1), the com- 
ponents of the radius and velocity are determined. For JUM?>1, the radius 
and velocity vectors are initialized by plumbline-state components stored in 
the input array VIV. After initialization of several variables, trigger 
values are set for the integration package for print, thrust events, tilt-over 
time, discrete altitudes, maximum dynamic pressure, Mach number equal to 1, 
and acceleration limits. Initialization of the trajectory is performed by 
first calling the subroutine ATHREV to determine the vehicle geometry, thrust 
levels, propellant flowrate, aerodynamics, integration steps, thrust event 
triggers, etc. ADER1 and ADER respectively, are called to determine the initial 
values of derivatives. The subroutine APRTN is then called with the header 
title "liftoff 1 *. 


The integration package is called and the trajectory Is integrated to 
termination, A check for integration error (KERR^O) is made. If no error has 
been made, the program flow is returned to the calling routine, otherwise, an 
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error message "INTEGRATION ERROR DUMPED" is printed. 


Functional Flow Diagram of AF0RUN 
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6.3.6 AGE0 

Subroutine Identification 


• Title 
AGE0 

• Calling sequence 
CALL AGE0 


Function 

Evaluates oblate gravitational properties 
Functional Flow 

This subroutine is used for calculation of terms necessary for determining 
the gravitational accelerations when the oblate model is required. 

Functional Flo w diagram of AGE0 
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6.3.7 AINIT 

Subroutine Identification 

• Title 
AINIT 

• Calling Sequence 
CALL AINIT 


s 


V 

« 


Function 

AINIT reads input data, and perforins initialization. 

Functional Flow 

Subroutine AINIT initializes variables and calls for NAMELIST input. 

After calling the input, specific values are converted to MKS system and most 
angles converted to radians. The number of optimized values and the number 
of constraints are checked and counted to insure being within the current 
program limits (ten parameters and ten constraints). Initial (launch or jump- 
start) values of the state variables and environmental terms are determined. 
The NAMELIST input is then printed for reference. 

Functional Flow Diagram of AINIT 
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6.3.8 AKALT 

Subroutine Identification 

• Title 

AKALT (Entry point in ATILT) 

• Calling sequence 
CALL AKALT 


Function 

Trigger routine used for 10 km and 14 km altitude printout. 

Functional Flow 

This subroutine is called by the integration package when the vehicle 
altitude reaches 10 and 14 kilometers. The subroutine APRTN is called to 
print the state variables with special printouts of "10 KMS." or "14 KMS." 

Functional Flow Diagram of AKALT 
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6.3.9 AMACH 

Subroutine Identification 


• Title 

AMACH (Entry point in ATILT) 

• Calling sequence 
CALL AMACH 


Function 

Trigger routine for MACH-1 printout. 

Functional Flow 

This subroutine is called when the Mach number is equal to 1. The 
subroutine APRTN is called to pric. the state variables with special printout 
of " M ACH ONE". 

Functional Flow Diagram of AM A CH 
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6.3. >0 AMULG 

Subrout ine Identification 


• Title 
AMULG 

• Calling sequence 

C: r L AMULG (L,M,N,X,XT,Y1, Y1T,Y2,Y2T), where 
L - number of dependent variables (1 or 2) 

M - previous value used in table location 

N - number of points in table 

X - independent variable 

XT - independent variable table 

Y1 - first dependent variable 

Y1T - first dependent variable table 

Y2 - second dependent variable 

Y2T - second dependent variable table 

Function 

Linear interpolation scheme 

Functional Flow 

(see appendix B) 

Functional Flow Diagram of AMULG 
(not required) 


6-20 






6*3.11 ANEWCH AND NNEWCH 

Subroutine Identification 

• Title 

ANEWCH; note that NNEWCH is an entry point in ANEWCH. 

• Calling sequence 

CALL ANEWCH or CALL NNEWCH 


Function 

Evaluates the parameter update equations to find the required changes 
in the controlling parameters. 

Functional Flow 

This subroutine uses the influence coefficients, constraint errors and 
the payoff values to determine the amount of change required in each controlling 
parameter to meet the constraints and optimize the payoff. 

Initially, the matrix of influence coefficients is adjusted such that 
those required in the optimization are relocated to the upper left portion of 
that total matrix. If two QY searches (see Section III) have not been 
initiated (indicated by N0SRCH < 2) a weighting matrix is determined which 
allows all variables to have similar weights; otherwise, the previous weighting 
matrix is used. The 1^ matrix is calculated as the product of the transpose 
of the constraint influence coefficient matrix, the weighting matrix, and the 
constraint influence coefficient matrix. The I fJ vector is calculated as the 
product of the first two matrices used to form 1^, and the payoff influence 
coefficient vector. 

The inverse of the I matrix is calculated by calling MATINV and the 
result is m?de symmetric by averaging across the leading diagonal. The product 
I ^ 1^ is found and the program flow proceeds to find the amount of con- 

trolling parameter correction desired to optimize the payoff and/or the amount 
desired to null constraint errors. The controlling parameter changes, dB, are 
then computed (see Section III) and added to the appropriate parameters. A 
test for convergence is made. If the desired convergence indicator tolerances 
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are satisfied, the required change In each parameter Is small, and LSB Is not 
equal to 1, the logical variable "BETC0N" Is set TRUE and the integer N0M1 is 
set equal to 1. If these conditions are not met, "BETC0N" is set FALSE 
and control Is returned to the calling routine. 


The entry point NNEWCH is called when a one-dimensional search is 
required, which does not require recalculation of the matrix products. 


Functional Flow Diagram of ANEWCH and NNEWCH 
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6.3.12 APRTN 

Subroutine Identification 


• Title 
APRTN 

• Calling Sequence 
CALL APRTN (NN), where 

NN - specifies the extra identification printout. 


Function 

Performs output editing, prints data, and creates tape file used for 
output tables and plots. 

Functional Flow 

Subroutine APRTN is a print subroutine called by subroutines ATHREV, 

ATILT, AMACH, AQMAX, AXPRT, GLINT, AF0RUN, and AKALT. If the subroutine is 
called during an influence coefficient determination run (LSB*1) , an immediate 
return is made; otherwise, the program checks the value of LPRINT (see Table 6 
for an explanation of the LPRINT options) to determine the amount and frequency 
of printed information desired. The subroutine calculates the variables 
desired in the printout which have not been defined through common block 
storage and converts those which have to the desired units of output. The 
print is made with a label determined by the variable NN. A check is made 
of the value of NTABLE to determine if tape output is desired. If not, the 
program returns to the calling subroutine. If tape output is desired, the 
pertinent data is stored on tape unit 9. 
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6.3.13 AQMAX 

Subroutine Identification 

• Title 

AQMAX (Entry point in ATILT) 

• Calling sequence 
CALL AQMAX 


Function 

Trigger routine for maximum dynamic pressure printout. 

Functional Flow 

Subroutine AQMAX is a print trigger routine ca'.led by the integration 
package when the time derivative of the dynamic pressure is zero, i.e., at 
the point where maximum dynamic pressure occurs. The routine turns the 
trigger off, and prints the state variable with a special printout of "Q 
MAXIMUM”. 


Functional Flow Diagram of AQMAX 
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6.3.14 ASIMP 

Subroutine Identification 


• Title 
AS IMP 

• Calling sequence 
CALL AS IMP 


Function 

Determines the flyback range using empirical data interpolation. 

Functional Flow 

The subroutine ASIMP is used to determine the flyback range based on the 
altitude, relative velocity, and relative flight-path angle at booster staging. 
A trivariant interpolation* is used to determine the flyback range based on 
empirical data. 

Functional Flow Diagram of ASIMP 



*See Appendix B for further information concerning 


interpolation routines. 
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6.3.15 ATHREV 
Subroutine Identification 

• Title 
ATHREV 

• Calling Sequence 
CALL ATHREV 

Function 

Thrust event control routine, sets up vehicle geometry, aerodynamics, 
number of engines, thrust levels, flowrates, and trigger flags for thrust 
event cutoff. 

Functional Flow 

Subroutine ATHREV is called at the beginning of each thrust event either 
by AFfifRUN before the integration starts or by the integration package during 
the trajectory run. This subroutine sets up the vehicle geometry, aerodynamics, 
number of engines, thrust levels, flowrates, etc., for the ensuing thrust 
event. ATHREV also stops the Integration at the end of the last thrust event 
by calling the subroutine RTMRK. 

ATHREV first determines which thrust event is baing initiated from the 
value of ITHR. If the first event is being Initiated (liftoff) ITHR will be 
equal to one and the first stage center-of-gravity, engine gimbal position, and 
aerodynamic data are entered into the tables used in the calculations, and the 
first stage pitch and yaw attitude polynomials are determined. If ITHR is 
greater than one, the value of LSTGE (ITHR) is checked to see if it differs 
from toe previous thrust event (ITHR-1), indicating that staging has occurred 
at the end of the last event. If staging has occurred, the second stage center- 
of-gravity, engine gimbal, position, and aerodynamic data are entered into the 
working tables, the flyback fuel required is calculated if desired,* and the 
second stage pitch and yaw attitude polynomials are determined. Also, the 
state aL staging 1 m Mlnred M the current trajectory Integration l*t not |>,n! <>l 
the influence coefficient calculation or the optimization step search. 

A See Appendix B for further information concerning flyback calculations . 
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Printout is required if the value of ITHR is not equal to (JUMP-1) . A test 
is then made to determine if the previous thrust event was the last event re- 
quired (by the variable NVNT) . If che final thrust event has been completed, a 
check is made to determine if the routine is being executed during an influence 
coefficient trajectory (LSB=1). If not, the subroutine calls APRTN for a final 
printout using the header "INJECTION" and stores the state derivatives. 

Whether or not the LSB variable is equal to 1, the final cutoff weight is cal- 
culated and the subroutines ADER1, ADER, and APRTN are called for a final pass 
through the equations of motion and a final printout. A call to the routine 
RTMRK terminates the integration package to return to the subroutine AF0RUN. 

If the final thrust event has not been completed, selection of cutoff 
triggers is made by the input variable MSWCH (see input section). Integration 
step size, print step size, exit area, propellant flowrate and thrust are 
calculated for the thrust event and a test is made to see if ITHR is equal to 
1. If ITHR*]. the program flow is returned to the calling routine, AF0RUN. 

If ITHR is greater than 1, the value of the previous thrust event, JT, is 
checked to see if an intermediate equality constraint is required at the end 
of the last event, noted by NVRST. If JT=NVRST, the state variables are 
stored at that time. The time derivative of dynamic pressure is then checked 
by calling ADERl and ADER to see if the ne* thrust event will result in a 
change in sign of the derivative. If there is a sign change, the subroutine 
AQMAX is called to mark maximum dynamic pressure. If the variable LSB»1 
(indicating influence coefficient trajectory), the program flow is returned 
to the calling routine. If LSB^O, the calculations for flight performance 
reserves are made if required and if the thrust event (ITHR) is equal to the 
input variable IPR. 

If the number of the previous thrust event is equal to NVRST, the state and 
derivatives of the state are stored and APRTN is called with a header "INJECTION". 
Whether the number of the previous thrust event was equal to NVRST or not, the 
subroutines ADERl and ADER are called to update the equations of motion. The 
APRTN subroutine is called with a header "THRUST EVENT". 
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Functional Flow Diagram of ATHREV 

















6.3.16 ATILT 

Subroutine Identification 


• Title 
ATILT 

• Calling Sequence 
CALL ATILT 


Function 

ATILT is a trigger routine for end of vertical rise (begin tilt- ever or 
attitude control) and end of tilt over (if used). the subroutines AKALT, 

AMACH, AQMAX, AXPRT, and GLIMT are contained within this subroutine as entry 
points . 

Functional Flow 

Subroutine ATILT is a print and control trigger routine that is called 
by the integration package at the end of the liftoff phase of the trajectory 
(specified by the input variable TLIFT) and again at the end cf the tilt-over 
phase if ang.le-of-attack attitude control is used. 

The subroutine first sets KS1=2 denoting end of vertical rise and turns 
off the ATILT trigger. A test is then made to determine if a programmed tilt- 
over is to be used (in conjunction with angle-of-attack attitude control). 

If it is, the ATILT trigger is turned back on and the trigger time is s 
the time at the end of the tilt-over, TTILT. If the call to ATILT was 
the end of tilt-over, the ATILT trigger is turned off, APRTN is called with 
the header "END TILT" and flow returns to the integration package. If a 
programmed tilt is not used, or if it is and the call to ATILT occurred at the 
end of the vertical rise, APRTN is called with uie header "BEGIN TILT", the 
maximum dynamic pressure trigger is turned on, the Xp and Xy attitude polynomials 
are determined by calling subroutine FIND, and flow returns to the integration 
package. If the call to ATILT has occurred at the end of a programmed 
tiltover, the ATILT t-igge is turned off, subroutine APRTN is called with 
the header "END TILT", and flow returns to the integration package. 
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Functional Flow Diagram for ATILT 


■*fc. 
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6.3.17 AXPRT 

Subroutine Identification 

• Title 

AXPRT (entry routine in ATILT) 

• Calling Sequence 
CALL AXPRT 


Function 

Trigger routine for normal printout. 

Functional Flow 

Subroutine AXPRT provides a call to the print routine APRTN during the 

trajectory integration indicated by the input variable LPRINT. The sub- 
routine first calls APRTN with a blank header and updates the value oi the 
next print time required. Using the function AINT the decimal portion o*. time 
is truncated such that print time will occur at whole number times. 

Functional Flow Diagram of AXPRT 
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6.3*18 BADLX 

Subroutine Identification 


• Title 
BADLX 

• Calling sequence 
CALL BADLX 


Function 

BADLX determines the influence coefficients (partial derivatives) used in 
the computation of the changes in the controlling parameters. 

Functional Flow 

Subroutine BADLX is used to determine the influence coefficients (the 
partial derivatives of the payoff and constraints with respect to the con- 
tolling parameters) which are used in the calculation of the changes in the 
controlling parameters. The subroutine first sets the flag LSB*1, which is 
used in other routines to identify the type of Integration run being made. 

The increment amounts for the parameters are determined and used according to 
the parameters required by the input KDB. A call to the subroutines AF0RUN 
and AF0RND is executed for each independent change of the controlling param- 
eters. The influence coefficients are then stored in the array XLAMB for 
use in ANEWCH. If the variable NPSTST is equal to 1, the printout "FORWARD 
DIFFERENCES” is printed, indicating that only one change will be made in the 
parameters. If the variable NPSTST is not equal to 1, the printout "CENTRAL 
DIFFERENCES" is printed which indicates that a plus and minus increment will 
be made in each parameter and the associated influence coefficients taken 
as the average of the plus and minus values. 

An exception to the above procedure for obtaining the partial derivatives 
occurs when final bum time is an optimized parameter and payload is not 
a payoff or constraint. In this case, the subroutine ACST0P provides BADLX 
a set of analytical partial derivatives based on the partial of the constraint 
with respect to the state variables. BADLX then determines the influence 
coefficients of the constraints with respect to bum time without integration 
of a trajectory simulation. 
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Functional Flow Diagram of BADLX 





6-37 


J 

* 

















6.3.19 B0PPLT 

Subroutine Identification 

• Title 
B0PPLT 

• Calling sequence 
CALL B0PPLT 


Function 

B0PPLT generates output plots on CALCOMP plotter 
Function of Flow 

The subroutine B0PPLT will be called from MAIN after completion of a 
converged trajectory by setting the input variable NTABLE equal to 2 or 3. 

This subroutine uses special CALCOMP library functions during execution. The 
programmer should verify that these routines are on the system library before 
execution of B0PPLT. 

The subroutine starts by reading a three-card, 27-word, field identifica- 
tion message used for communication with the plot opere*' . The identification 
message is printed and a call is made to the CALID routine which plots an 
identification print and begins the plotting process. The PLOT subroutine is 
called which sets up the location o £ the origin in terms of location on the 
plot. Plotting information is supplied by an input card which determines the 
variables required for plotting and the desired units for output. If the word 
ENDPLT is input, the CALEND subroutine is called which draws a final identi- 
fication block and terminates the plotting routine. 

For each input card the variables required and the units desired are 
identified and the correct variables are taken from the trajectory tape pro- 
vided by APRTN. These variables are read from the tape and placed into 
storage arrays. The SCALE subroutine is called which determines the scaling 
required for the plotting process. 
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The AXIS subroutine is called which labels the plot and draws the axes. 
The subroutine LINE is then called to draw the generated curve. 

If labeling is required* a label is read from cards and the subroutine 
SYMB0L is called. Each plot is concluded by calling the PL0T subroutine to 
establish a new origin. Program flow returns to read another plot variable 
input card. 


Functional Flow Diagram of B0PPLT 














6.3.20 B0PTPL 

Subroutine Identification 


• Title 
J0PTBL 

• Calling Sequence 
CALL B0PTBL 


Function 

B0PTBL generates output table summary. 

Functional Flow 

The subroutine B0PTBL will be called from MAIN after completion of a 
converged trajectory by setting the input variable NTABLE equal to 1 or 3. 

The subroutine uses a special NAMELIST call, therefore additional input is 
required when using this subroutine. 

The general flow of the subroutine is to read data from a trajectory 
tape provided by the APKTN subroutine and place this data in the fourteen 
output tables. These tables (see subsection 9.3) are comprised of two sets of 
seven tables each which provide output in the MELS system and the English 
system. 
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6.3.21 • CHIP0L 

Subroutine Identification 

• Title 
CHIP0L 

• Calling sequence 

CALL CHIP0L (N1,N2,T,Y1,A1,Y2,A2) where 

N1 - order of A1 polynomial 

N2 - order of A2 polynomial 

T - independent variable 

Y1 - value of polynomial Al evaluated at T 

A1 - polynomial coefficients 

Y2 - value of polynomial A2 evaluated at T 

A2 - polynomial coefficients* 

Function 

Evaluate the pitch and yaw attitude polynomials. 

Functional Flow 

(not required) 

Functional Flow Diagram of CHIP0L 
(not required) 
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6.3.22 FIND 

Subroutine Identification 


• Title 
FIND 

• Calling Sequence 

CALL FIND (CHI,A,DT,N0RDER) where 
CHI - attitude table 
A - output polynomials 

DT - time table based on beginning of stage time 
N0RDER - order of the polynomial desired. 

Function 

FIND determines attitude polynomials based on tabulated attitude time 
history information. 

Functional Flow 

(not required) 

Functional Flow Diagram of FIND 
(not required) 
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6.3.23 GLIMT 

Subroutine Identification 

• Title 

GLIMT (Entry point in ATILT) 

• Calling Sequence 
CALL GLIMT 


Function 

Trigger routine for acceleration limit printout. 

Functional Flow 

The subroutine GLIMT is called from the integration package when the 
axial acceleration limit has been reached. 

The acceleration trigger is turned off and the APRTN subroutine is called 
with the header "BEGIN GLIMIT”. If discrete throttling is required (MSWCH<0) 
the subroutine ATHREV is called; otherwise, the flow is returned to the 
calling routine. 

Functional Flow Diagram of GLIMT 
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6.3.24 MATINV 

Subroutine Identification 

• Title 
MATINV 

• Calling Sequence 

CALL MATINV (B,N) where 

B - matrix to be inverted (the inverted matrix is stored in B also) 
N - number of rows of matrix B. 


Function 

MATINV performs matrix inversion, calculating the inverse in double 
precision arithmetic, although the entry and exit matrices are truncated to single 
precision. 

Functional Flow 

(not required) 

Functional Flow Diagram of MATINV 
(not required) 
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6.3.25 PRA63 

Subroutine Identification 


• Title 
PRA63 

• Calling Sequence 

CALL PRA63 (PR.ERR0R) where 

PR - storage array for atmospheric parameters 

ERR.0R - error flag (not used). 


Function 

PRA63 evaluates the atmospheric parameters as well as their partial 
derivatives w.r.t. altitude using the spline interpolation method of the 
subroutine SPLINE. 

Functional Flow 

Interpolation method used in this routine is defined in appendix B. 

Functional Flow Diagram of PRA63 
(not required) 
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6.3 ?6 SPLINE, SPLIN2 
Subroutine Identification 


• Title 
SPLINE, SPLIN2 

• Calling sequence 

CALL SPLINE (L,M,N,X,XT,Y,TY,YD,K,LAST) 

CALL SPLIN2 (L,M,N,X,XT,Y,TY,YD,K,LAST) where 
L - number of dependent variables (stored in comw>ri storage) 

M - location index 

N - number of points in table 

X - independent variable 

XT - independent variable table 

Y - dependent variables 

YT - dependent variables tables 

YD - first derivative of the dependent variable w.r.t. the independent 
variable 

K - flag to denote calculation of partial derivatives 

K=l, do not calculate partials; K*2 calculate partials 

LAST - stored value of M. 


Functi. i 

SPLINE determines the aerodynamic coefficients using a SPLINE interpolation 
method. 

SPLIN2 is identical to SPLINE in logic, however, SPLIN2 is used for 
interpolation of thrust and weight loss for the SRM engines. 


Functional Flow 

(see appendix B) 

Functional Flow diagram of SPLINE 
(not required) 
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6.3.27 TEST 

Subroutine Identification 


• Title 
TEST 

• Calling Sequence 

CALL TEST (X,XT,NX,I,R) where 

X - Independent variable 

XT - independent variable table 

NX - number of points in XT table 

* - location index 

R - dependent variable slope 


Function 

TEST embodies the linear schenu used in ASIMP for determining slope of 
the dependent variable w.r.t. the independent variable. 

Functional Flow 

(not required) 

Functional Flow Diagram of TEST 
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4.12B TRIVLP 
Subroutine Identification 


• Title 
TRIVLP 

• Calling Sequence 

CALL TRIVLP (G,GT,NG,V,VT,NV,A,AT,NA,<|>T,<J.) where 

G - value of relative flight-path angle 

GT - relative flight-path angle table 

VG - number of points in GT table 

V - value of relative velocity 

VT - relative velocity table 

NV - number of points in VT table 

A - value of altitude 

AT - altitude table 

HA - number of points in AT table 

4>T - flyback range table 

$ - flyback range (output). 

Function 

Trivariant lookup of flyback range (for use in ASIMP) . 

Functional Flow 

The subroutine TRIVLP is called by the subroutine ASIMP to determine the 
flyback range as a function of relative flight-path angle, relative velocity, 
and altitude. 

ASIMP calls the subroutine TEST to determine the linear slopes of the 
altitude, relative velocity, and relative flight-path angle w.r.t. the 
flyback range. The flyback range is then calculated using these slopes. 

Functional Flow Diagram of TRIVLP 


(See TEST) 



PART III 
USER S MANUAL 



Section VII 

PROGRAM LISTING 

Consult MSFC, S&E t AERO-GT for Program Listing 
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Section VIII 

VARIABLE NAME CROSS REFERENCE OF MAJOR SUBROUTINES 


Consult MSFC, S&E, AERO-GT for variable name cross reference of major 
subroutines. 


8—1 



Section IX 

INPUT /OUTPUT 


The usefulness of any computer program as an analysis tool is dependent 
largely upon how well the user understands the operation of the program and 
the input and output involved in its use. This section has therefore bet., 
designed to provide the user a thorough understanding of the basic mechanics 
of tVe program, the input required, and the output available. The ^uwing 
paragraphs describe in detail; 1) how to formulate the trajectory problem for 
presentation to RAGMOP, 2) how to select the various options available, 3) the 
physical input data required. 4) the various types of output produced by the 
program, and 5) some suggestions concerning technique which are helpful in 
making the most efficient use of the program 

RAGMOP users familiar with the ROBOT computer program will note the 
similarity in the input of the two programs. The input philosophy of RAGMOP 
is intended to conicide with that of ROBOT whenever possible. This was done 
(due to the wide use of ROBOT) as a convenience to potential users of RAGMOP. 
Note, however, that beyond this similarity the two programs are vastly dif- 
ferent. 

9.1 PROBLEM FORMULATION 

This section describes how to formulate the ascent trajectory for 
presentation to RAGMOP. The terminology used in describing various events 
is defined, and information pertinent to the selection of various program 
options is presented. 

9.1.1 Thrust Profile 

The RAGMOP trajectory is described as a succession of a number of thrust 
events, one or more of which may comprise a complete stage. A thrust event 
is a period of time characterized by a continuous thrust profile (constant, 
zero, or varying in a continuous fashion) (Figure 9-1) . Each thrust event is 
initiated at the termination of the previous event except, of course, for 
the first event which initiates the trajectory. A thrust event is terminated 
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Figure 9-1. THRUST EVENT PROFILE 


on one of four criteria: (1) time (the duration of the thrust event, not the 

absolute final time), (2) liquid fuel depletion (fuel is input in the first 
thrust event of a liquid bum and any fuel remaining at the end of a given event 
is carried forward to the next event of the same stage, if any), (3) g-limit, 
or (4) relative velocity (which, of course, should be used only for one thrust 
event). The data used to determine the vehicle geometry, propulsion, and aero- 
dynamics is described either by thrust event or by stage. Data which is input 
by stage consists of the aerodynamic coefficient, base drag, and center-of- 
gravity tables, the aerodynamic reference length, the order of the Xp and x Y 
attitude polynomials, number of thrust events comprising the stage, engine 
gltnbal point locations, and the location of the aerodynamic data moment ref- 
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ereace point. Data which is input by thrust event consists of the number of 
engines, vacuum thrust per liquid engine, exit area per SRM or liquid engine, 
mass flowrate per liquid engine, specific impulse for liquid engine when con- 
tinuous throttling is used, liquid fuel per thrust event (actually input in the 
first thrust event of the liquid stage which begins the liquid burn), accelera- 
tion (g) limit, liquid thrust event cut-off flag, print and integration step sizes, 
aerodynamic reference area, thrust event duration, and jettison weight (which 
is dropped at the end of the thrust event) . 

Data for the SRM engines consisting of sea level thrust and weight loss 
are input as tabular functions of time from liftoff. 

9 1.2 Control Program 

In formulating the trajectory problem, the user must decide what kind of 
attitude control Is required* If the optimized polynomial form is used for the 
Xp program, the degree of the polynomial desired for each stage must be specified. 

Also, the decision must be made as to whether or not the yaw program is to be 

optimized, and if so, the degree of the yaw polynomial in each (or either) stage 
must be selected. If angle-of-attack control is to be used for the first stage, 
the optimized tilt-over maneuver must be used. This requires selection of the 

polynomial degree and the length of time desired for the pitchover. 

Note that a polynomial for Xp and Xy can b® used without being optimized 
or without each data point used in the curvefit being optimized. Suppose, for 
example, that the value of Xp at a certain time during the first stage was re- 
quired to be a fixed value, and that a cubic polynomial was desired for the Xp 
of the stage. The cubic polynomial could be used with only two optimized data 
points, with the first point and any one of the other three points (four points 
required to perform cubic curvefit) remaining fixed throughout the run. As an- 
other example, a Xy polynomial could be used without optimization of any of the 
data points. 

9.1.3 Efficient Use of RAGMOP 

While a great many combinations of parameters and constraints exist within 
RAGMOP, careful selection of the optimized parameters and the constraints en- 
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forced can streamline the program greatly with little or no effect on the s^iu- 
tion obtained. Several comments and suggestions are presented below which pro- 
vide information pertinent to using RAGMOP most efficiently. 


9. 1.3.1. Body axis Coordinate System . In order to agree with the present 
Apollo Standard coordinate system, all vehicle input data uses the coordinate 
system shown in Figure 9-2, with the longitudinal axis as the x-axis, the lat- 
eral axis as the y-axis (out the right wing), and the vertical axis as the z- 
axis (down). At liftoff, the vehicle will be oriented with the z-axis pointed 
in the direction of the launch azimuth (downrange) . 


9.1. 3.2 Constrained qa and/or qB . Although RAGMOP has the capability of en- 
forcing the inequality qa and qB constraints, experience has shown that often 
these constraints will not be violated when a q^^ constraint is simultaneously 
enforced. Since the use c ' these constraints presents discontinuities in the 
trajectory angle-of-attack and/or sideslip angle derivatives, the solution using 
the constraints is generally more difficult to obtain and hence more computer 
time will be required. If some doubt exists as to whether or not the qa and/or 
qB limits will be exceeded, a solution can first be obtained without enforcing 
the constraints. If the maximum values are exceeded, the solution may be used 
as the initial "guess’* (nominal) trajectory for a run with the constraints enforced. 


9.1.3. 3 Constrained q^^. The q^ ax constraint in RAGMOP is enforced as an 
equality constraint. Therefore, as with the qa and qB constraints mentioned 
above, the constraint should not be enforced unless the user is sure that the 
maximum allowable value will be exceeded. Since the q^^ constraint is an equal- 
ity constraint, the optimized parameters will be adjusted until the constraint is 
met. Thus, all parameters (launch weight, aximuth, burn times, and attitude con- 
trol) are affected. Experience has shown that creating a thrust event prior to 
^raax wliere t ^ irust is reduced by a discrete amount (at an optimized time) will 
aid in satisfying the constraint with a minimum of reshaping of the trajectory. 
This combination of throttling and reshaping has been seen to be much more ef- 
ficent than reshaping alone. 
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9. 1.3.4. Payload Option . The payload calculation (see Appendix D) has been 
designed especially for space shuttle type vehicles and is therefore preferable 
to the FPR option for these vehicles. 


9. 1.3.5. Attitude Control Program . The vehicle attitude control polynomials 
are stage-dependent and therefore are carried over thrust events of a given 
stage. The order of the polynomial used to describe the attitude histories 
should be chosen with the following points in mind: 

• The higher the order of the polynomials, the more optimum the trajec- 
tory can be. 

• The higher the order of the polynomials, the more parameters are re- 
quired and therefore the longer the run will take to converge. There- 
fore, the lowest order polynomial which yields an acceptable solution 
should be used when rapid solutions are required. 

• An optimized launch azimuth or yaw control program is generally not re- 
quired unless a specific orbital inclination is a constraint, or when 
the coordinated turn option is used. 

• The second stage Xp program is normally initiated at the last value of 

Xp from the first stage (continuous Xp) • The use of a discontinuous 
Xp is optional (KDB(19)=1) and selection of discontinuous Xp can be 
based largely on two points: (1) the discontinuous x p is in general 

more optimum and (2) the use of Xp discontinuity at staging may not be 
realistic if the dynamic pressure is above a certain level (i.e., dif- 
ficulties in executing the Xp discontinuity may arise in actual flight). 
If the second stage is flown essentially in a vacuum, (i.e. no aero- 
dynamic coefficients input) there is no need for more than a linear x p 
profile since, even with a quadratic form available, the Xp program 
produced the RAGMOP will be linear (due to the nature of the vacuum 
flight problem). 

In addition to the above points concerning selection of the order of the 
attitude polynomials, one should also note that proper spacing of the points in 
the first and the second stage time tables can help increase the efficiency of the 
program. The input tables TTBL (first stage) and TOBL (second stage) contains the 
time points at which the angles in the input tables CPTBL, CYTBL (first stage), 
CP0TBL, and CY0TBL (second stage) apply. The Xp and Xy values in the latter set 
of tables will be varied at each parameter update step, and curvefits for the 
appropriate polynomials performed, in light of this, it is easily seen that an 
even distribution of time points over the duration of each stage will provide 
the most stable operating conditions for the program. (Consider, for example, 
the difficulty of fitting a curve to points grouped very closely together, and 


the sensitivity of the curve to very slight variations in the points when they 
are so grouped) . 

9. 1.3. 6 Run Time . It should be noted that RAGMOP, while appearing similar in 
some ways to ROBOT, may at times require considerably more computer time to 
obtain a converged trajectory solution. This is due to the more complex physi- 
cal model (atmosphere in both stages, optimized from launch to orbit, lifting 
trajectory, etc.) and to the necessarily different optimization scheme. If a 
particular run is returned and has not obtained a converged solution, the user 
may take advantage of the work performed during the first run ty merely updating 
the values of the optimized parameters and resubmitting the run. Generally, 
the parameter set obtained during a run will meet the required constraints, but 
if the run is terminated due to maximum time, will not be optimum. Thus, a 
run which is not yet converged may still provide the user with meaningful in- 
formation, unlike other techniques wherein intermediate results have no signifi- 
cance. Trajectory solutions with RAGMOP have been obtained in from two to ten 
minutes of computer time starting from very poor nominal trajectories. Good 
initial guess (nominal) trajectories will increase speed of convergence. 




Figure 9-2. INPUT BODY AXIS SYSTEM 
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9 2 ITEMIZED DESCRIPTION OP NAMELIST INPUT 


NAME 

DIMENSIONS 

UNITS 

EXPLANATION 

AE 

(15) 

M 2 

Exit area per engine for each thrust event. 
Preset to zero. 

ALAT 


DEG 

Launch site latitude. Preset to 28. 531885. 

ALONG 


DEG 

Launch site longitude. Preset to 80.5649528. 

ALTBAS 

(25,2) 

M 

Altitudes at which the base pressure ditter- 
entials in BAXIAL apply. Up to 25 altitudes 
may be used. The second index denotes the 
combined booster-orbiter (1) or the orbiter 
alone (2). Preset to zero. 

ALTLS 


M 

Altitude of the launch site above the spher- 
oid. Preset to zero. 

ALTTBL 

(25) 

M 

Altitude values at which the wind speeds and 
directions in WTBL and AZWTBL apply. 

AYL 



Used for error check in integration. The 
preset value should be used. (No input re- 
quired). Preset to 0.002. 

AZ 


DEG 

Launch azimuth. Preset to 90. 

AZWTBL 

(25) 

DEG 

Azimuth angles denoting the directions that 
applv to the wind speeds and altitudes in 
WTBL and ALTTBL. Denotes the direction 
toward which the wind is blowing. 

BAXIAL 

(25,2) 

LBF 

Base axial force. 

BL0O 



S/S) b° oster * Used in booster flyback 

fuel computation. 

CAALP 

(25,2) 

/DEG 

Partial of axial force coefficient with res- 
pect to angle-of-attack. Indices as in CA0. 
Preset to zero. 

CA0 

(25,2) 


Zero lift axial force coefficient correspond- 


ing to Mach numbers in PNM. Second index 
denotes stage. Preset to zero. 
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NAME 

DIMENSIONS 

UNITS 

EXP LAW AT ION 

CASK 



The number ol the rase >e Lng run. For output 
use oniy. Preset to 1. 

Cl SI’ 

(i3) 

SEC 

Vacuum L of liquid engines used during 

sp 

continuous throttl ing. 

C ) 



First cue! 1 icient Ln the Fischer ellipsoid 
gravitational expansion. The preset value 
should be used. (No input required) 

C LB ETA 

(25,2) 

/DEG 

Partial of rolling moment eoetlicient with 
respect to sideslip angle. indices as in 
CA0. Preset to zero. 

CMALP 

(25,2) 

/DEG 

Partial of pitching moment coefficient with 
respect to angle-of-attack. Indices as in 
CA0. Preset to zero. 

CM0 

(25,2) 


Zero lifL pitching moment coeflieient. in- 
dices as Ln CA0. Preset to zero. 

CMUK 


3 2 

M /SEC 

Product of the universal gravitational con- 
stant and the mass of the earth. The pre- 
set value should be used. (No input required) 

CNALP 

(25,2) 

/DEG 

Partial of normal force coefficient with 
respect to angle-of-at tack. Indices as in 

CA0. Preset to zero. 

CNBETA 

(25,2) 

/DEG 

Partial of yawing moment coefficient with 
respect to sideslip angle. Indices as in 
CA0. Preset to zero. 

CN0 

(25,2) 


Zero angle of attack normal force coefficient. 
Indices as for CNALP. Preset to zero. 

C0RBWT 


LBM 

Constant orbiter weight used in payload 
calculation. 

C.P0TBL 

(20) 

DEG 

Values of chi-pitch (y ) to be used in form- 
ing the nominal trajectory for the second 
stage* Values correspond to time in T0BL(l-iO) 

CPTBL 

(20) 

DEG 

Same as CP0TBL but for first stage, and 
values correspond to times in TTBL(l-JO) 

CRSDWT 


LBM 

Constant residual weight used if SRKSI D = 0. 
Preset to zero. 





r 
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NAME 


DIMENSIONS 


UNITS 


EXPLANATION 


CTNKWT 


LBM 

Constant tank weight used if SCALE « 0. 
Preset to zer^ 

CWD0T 

(15) 

LBM/SEC 

Critical flow rate per engine for each thrust 
event . 

CYBETA 

(25,2) 

/DEG 

Partial of side force coefficient with res- 
pect to sidelsip angle. Indices as in CA0. 
Preset to zero. 

CY0TBL 

(20) 

DEG 

Same as CP0TBL but for yaw (correspond to 
T0BL( 11-20)) 

CYTBL 

(20) 

DEG 

Same as CPTBL but for yaw (correspond to 
TTBLU1-20)) 

DELVG 


M/SEC 

Delta velocity required for geometry reserves. 
Preset to zero. 

DELVP 


M/SEC 

Delta velocity required for performance 
reserves. Preset to zero. 

DJ 



Third coefficient in the Fischer ellipsoid 
gravitational expansion. The preset value 
should be used. (No input required). 

DTZ 


SEC 

Time from ground reference release (GRR) to 
liftoff. GRR is the point in the countdown 
at which the launch inertial coordinate sys- 
tem is established. Must be input for jump 
starts as well as for ground-launch trajec- 
tories. Preset to zero. 

DV0MS 


FT/ SEC 

Delta velocity required of orbital maneuvering 
system, used in payload calculation. 

END 

(10) 


Maximum allowable absolute values of the 
constraint errors. END(I) corresponds to 
KCDPH 1(1+1) , and then to KCDRES(I+N) where N 
is the number of constraints in KCDPHI. This 
set of tolerances must be met before a run is 
considered to be converged. 

E0PCT 



Additional fuel percentage needed for booster 
flyback with one engine out. Percentage of 
landed weight of booster. 

EU 



Upper error bound in integration. The preset 


value should be used. (No input required). 
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NAME 

DIMENSIONS 

UNITS 

EXl’LANTle.i 

FACT 



Denotes whether optimised lifting (FACT = -I), 
zero normal force (FACT = 1), zero angle-of- 
attack (FACT * 2) or Mach number dependent 
(FACT = 3) is to be flown. If FACT is 1,2, 
or 3, a programmed tilt over time should be 
used (See TTILT). Preset to -1. 

FLBS 

(15) 

LBF 

Liquid engine vacuum thrust per engine for 
each thrust event. Preset to zero. 

FPRFAC 



Factor uoed in calculation of delta velocity 
for performance reserves (FPR) 

Cl (fiw 

FPR - FPRFAC* [g Q V1SP in ~~ + g Q VlSPO In 

0L0W , 

0L0W J 

FUELIQ 

(15) 

LBM 

Liquid engine fuel* For parallel burn, all 
liquid ermine fuel should be put Into FUELIQ 
(1) and the indicator NPARBN set equal to 1, 
For series burn, the fuel is irput in the 
first thrust event of each stage. Must be 
input whenever tuel consumed is used as a 
cutoff criteria ftr any thrust event r * a 
stage. 

GAFPCT 



Additional fuel percentage required by booster 
for a go-around at the end of flyback. Per- 
centage of landed weight booster* 

GLIM 

(15) 

G's 

Maximum longitudinal acceleration, in g*s for 
each thrust event. Preset to 0.1E20. 

GZER0 


M/SEC 2 

Gravitational acceleration of earth at the 
equator. The preset value should be used. 
(No input required). 

H 



Second coefficient in the Fischer ellipsoid 
gravitational expansion. The preset value 
should be used. (No input required). 

HEAD 

90 Spaces 


90 Space Hollerith Field used for identifi- 
cation of each page of trajectory print. 

HEDWND 


ft/sec 

Headwind used in fxyback fuel calculation. 
Preset to zero. 

HMN 



Minimum step-size for integration. The pre- 
set value should be used. (No input required). 
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NAME 


DIMENSIONS 


UNITS 


EXPLANATION 


IPR 

J0RB 

JTHR (15) 

JUMP 

KBACK 

KCDPHI (10) 

KCDRES (6) 

KDB (30) 

KDT (10) 


Denotes the thrust event from which the per- 
formance reserves are taken. Must be in the 
last stage and no thrust evert may be opti- 
mized after this one. If IPR is greater than 
zero, WPM and CWD0T must be input. Preset 
to zero. 

Denotes whether spherical (J0RB * 1) or ob- 
late (J0RB=O) earth model is desired. Pre- 
set to zero. 

Denotes use of SRM thrust and delta weight 
tables for e^ch thrust event (JTHR * 1) or 
no SRM thrus or delta weight (JTHR = 0), 

Preset to zero. 

Denotes the number of the thrust event at 
which the start is to occur. A normal ground 
launch would start at the first thrust event 
(JUMPED. A jump start can start at any later 
thrust event. For a jump start the array 
VIV must be input. Preset to 1. 

Indicator to denote computation of booster 
flyback fuel desired (KBACK = 1) or not 
desired (KBACK * 0). Preset to zero. 

Terminal constraint and payoff codes. 

KCDPHI(l) denotes the payoff desired. 

KCDPHI (2-10) are the code numbers of the 
terminal constraints desired. See Table 9-1. 
Preset to KCDPHI = 1,2 ,3, 4. 

Intermediate constraint codes. See Table 9-1. 
Preset to zero. 

Parameter optimization indicator. Each loca- 
tion corresponds to a particular control pa- 
rameter. (See Table 9-2 of this section. 

KDB (1)=1 denotes that the corresponding 
parameter is to be optimized. KDB(I)=0 denotes 
no optimization of that parameter. Preset 
to zero. 

A companion matrix to KDB(l-lO). Denotes 
the number of the thrust event from that in 
KDB(I) which is to be altered to hold tank 
limits. If the first four thrust events are 
to be optimized (i.e. , their burn times) with 
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NAME 


DIMENSIONS 


UNITS 


EXPLANATION 


KIND 

KP 

KRDER 

KY 

LAST 

LPRINT 

M0MBAL 

M.SWCH 

NC00RD 


the tank limits held in thrust event 
5 for the first two events and thrust event 
8 for the second two, input KDB=1, 1,1,1, 
and KDT=4,3,5,4. Preset to zero. 

Denotes the type of integration desired. 
KIND=1,2,3 denotes: variable step-size 

Adams-Moulton, Runge-Kutta, or fixed step- 
size Adams , respectively. KIND=1 is usually 
preferable. Preset to one. 

(2) Denotes the order of the chi-pitch polynomi- 

al desired for each stage. Preset to zero. 

Order of differences in the integration 
package. The preset value should be used, 
(no input required). Preset to 3. 

(2) Same as KP but for yaw. Preset to zero. 

Denotes only one data pack (case) to be 
evaluated (LAST=1) , or more than one (LAST= 
number of cases). n reset to one. 

Print option indicator. See Table 9-4. 
Preset to zero. 

Moment balance option indicator. Preset 
to zero. (M0MBAL = 0, no moment balance; 
M0MBAL = 1, SRM’s thrust through specified 
point (SRXCGF,SRZCGF) , liquid engines bal- 
ance moments; M0MBAL « 2, SRM’s track c.g., 
liquid engines balance moment; M0MBAL = 3, 
all engines balance moments is collective 
two engine equivalent model.) 

(15) Thrust event cutoff indicator. Tells the 

program which cutoff triggers are to be 
turned on for each thrust event. See Table 
9-3. Preset to 4 for each thrust event. 

(15) Coordinated turn option indicator per thrust 

event, NC00RD * 0, no coordinated turn; 
NC00RD = +1, coordinated turn with positive 
angle of attack; NC00RD * -1, coordinated 
turn with negative angle of attack. 
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NAME 


DIMENSIONS 


UNITS 


EXPLANATION 


NMAX 

N0EVNT (5) 

XPARBN 

NTABLE 

NVRST 

NWIND 

0MEGA 

0MSISP 

PNM (25,2) 

PRINT (15) 


Maximum number of trajectory integrations al- 
lowed in constraint-zeroing portion of pro- 
gram* Preset to zero* 

Number of thrust events in each stage, e.g., 
if there are two stages, with five thrust 
events in the first stage and two in the sec- 
ond,, input N0EVNT=5,2. Preset to cero. 

Parallel burn indicator. Set NPARBN=1 if 
parallel burn is desired, and place all fuel 
for orbiter engines into FUELIQ(l) (if 
fuel cutoff is required). If parallel burn 
is not desired, fuel for liquid engines must 
be input in the first thrust event of each 
stage for which fuel cutoff is required. 

Denotes output of tables and/or plots at end 
of converged run. Preset to zero. (NTABLE 

* 0, no post-processing is required; NTABLE 

* 1, output tables only; NTABLE -2, output 
plots only, NTABLE =3, output tables and plots.) 

Denotes the number of the thrust event at the 
end of which the intermediate constraints are 
to be imposed* Must be zero if no interme- 
diate constraints are desired. Preset to zero. 

The number of points used in the wind tables 
WTBL, AXWTBL, ALTTBL. If NW1ND=0 notables 
are used. Preset to zero. 

RAD/SEC Angular rotational velocity of the earth. The 
preset value should be used. (No input re- 
quired) . 

SEC Specific impulse of the orbital maneuvering 

system. 

Mach numbers at which the values of the vari- 
ous aerodynamic coefficients and the alpha 
history in TALFP apply. Indices as in CA0. 

((toy number of points may be used for either 
stage.) Preset to zero* 

SEC Print time increment for each thrust event 

during a long trajectory print. Preset to 
10 seconds for each event. 
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NAME 


DIMENSIONS 


UNITS 


EXPLANATION 


PSIREQ (10) 


PSIRST (6) 


OALPM 


QBETAM 

RE 

S (15) 


Desired values of the end constraints. PSIRKQ 
(1) corresponds to KCDPUL (1+1 ) . See Table 9-1 
Preset to 7876. 4195, .01 , 6470762. (V,y» and 

R for 50 & 100 n mi orbit). 

Desired values of the intermediate constraints 
PSIRST(I) corresponds to KCDRES(l). Preset 
to zero. 

2 

LB DEG/FT Maximum allowed value of the product of dyna- 
mic pressure and angle-of-attack. If the xp 
history causes this value to be exceeded, al- 
pha is reduced to hold the limit, overriding 
the x polynomial. Preset to 1.E20. 

2 

LB DEG/FT Same as QALPM but for sideslip angle (beta) 
and Y • Preset to 1.E20. 

y 

M Radius of the earth. The preset value should 

be used. (No input required). 

2 

M Aerodynamic reference area for each thrust 

event. Preset to zero. 


SCALE 

SCRSD 


Scale factor used in computing tank weight 
from propellant weight in payload calculation. 

Residual scale factor (decimal fraction). 
Preset to zero. 


SFC 

SPRADB 

SPRAD0 

SPMAE 

SRDWTB (15) 

SRMFTB (15) 


LB/LB-HR Booster flyback thrust specific f uel consump- 
tion. 

Perturbation value of Xp points used to form 
first stage Xp polynomial. Preset to .1°. 
Used in influence coefficient calculation. 

Same as SPRADB but for second stage. Preset 
to .5°. 

SRM exit area per engine 

LBM Total SRM weight overboard corresponding to 

time points in SRMTTB. (Note: not delta 

weight per engine.) 

LBF SRM sea level thrust per engine corresponding 

to time points in SRMTTB. 
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NAME 


DIMENSIONS 


UNITS 


EXPLANATION 


SRMPRP 


LBM 

SRMTTB 

(15' 

SEC 


SRXCGF 


M 

SRXGP 

(5) 

M 

SRYGP 

(5) 

M 

SRZCGF 


M 

SRZGP 

(5) 

M 

STEP 

(15) 

SEC 


SYRADB 



SYRAD0 



TALFP 

(25,2) 

DEG 

TAUT 

(15) 

SEC 

TBDWT 

(15,2) 

LBM 


TUFT 


SEC 


Total SRM fuel available. Should be the same 
as the last data point in the SRDWTB array. 

Time points at which the thrust levels and 
delta weights in SRMFTB and SRDWTB apply 
(note: a thrust event should be created 

whenever the slope of the thrust-time curve 
is discontinuous*) 

Longitudinal coordinate of SRM aiming point 
for fixed SRM’s (M0MBAL=1) 

Longitudinal coordinate of the SRM gimbal 
points. 

Lateral coordinate of the SRM gimbal points. 

Vertical coordinate of SRM aiming point for 
fixed SRM’s. 

Vertical coordinate of the SRM gimbal points. 

Integration step-size increment for each 
thrust event when Runge-Kutta or fixed 
step-size Adams-Moulton integration is used. 
Preset to 1 second for the first thrust event 
and 8 seconds thereafter. 


Same as SPRADB but for yaw. 

Same as SPRAD0 but for yaw. 

Angle-of-attack values corresponding to Mach 
numbers of PNM, Preset to zero. 

Duration of each thrust event. If a thrust 
event time is being optimized, this value is 
used as the initial estimate. Preset to zero. 

Delta weight values corresponding to the 
c.g. locations in TXCG and TZCG. Up to 
15 points may be used for each stage. 

The second index denotes the stage. 

Time at end of liftoff phase and the begin- 
ning of the Xp and/or Xy control program. 
Referenced to TZER0. Preset to 8 seconds. 
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NAME 


DIMENSIONS 


UNITS 


EXPLANATION 


TNE 

T0BL 

TTBL 

TTILT 

TXCG 

TZCG 

TZER0 

VCRUS 

VISPB 

VISP0 

VIV 


(4,15) 


Number of liquid engines per thrust event 
(TNE(1 ,ITHR) ) , number of SRM engines per 
thrust event (TNE(2,ITHR)) , and SRM pointing 
angles for fixed SRM's (TNE(3&4 ,1) ) • The 
latter two values are calculated internally 
and are not input requirements. 

(20) 

SEC 

Times, from staging, at which theXp and 
values in CP0TBL and CY0TBL apply. (T0BL 
(1-10), refer to TOBL (11-20) refer to 

X v > 

(20) 

SEC 

y 

Same as T0BL but for first stage and refer- 
enced from TLIFT. 


SEC 

Time at end of programmed tilt-over from 
which zero-lift, zero angle-of-attack, or 
Mach number dependent alpha profile is to 
be flown. The programmed tilt should be 
accomplished using a linear Xp (KP*1) for the 
first stage. One value of CPTBL and TTBL 
should be input, and the tilt may or may 
not be optimized, according to the coding 
in KDB. 

(15,2) 

M 

Longitudinal coordinate of the center of 
gravity corresponding to delta weights in 
TBDWT. See Figure 9-2. 

(15,2) 

M 

Vertical coordinate of the center of gravi- 
ty corresponding to delta weight' in TBDWT. 
See Figure 9-2. 


SEC 

Time at the start of the trajectory. For a 
normal launch TZER0 is the time at lift-off. 
For a jump start TZER0 is the time at the 
initiation of the trajectory. 


FT/SEC 

Booster flyback cruise velocity desired. 


SEC 

Vacuum 1 of booster. Used in payload cal- 

sp 

culation. 


SEC 

Vacuum I of orbiter. Used in payload cal- 

sp 

culation. 

(15) 


Initial state for a jump start. If VIV(7)=0, 
the plumbline state vector, W,U,V,X,Y,Z 
(£,£,?,Z,X,Y in the Apollo 13) must be input 
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NAME 


DIMENSIONS 


UNITS 


EXPLANATION 


into VIV(l-6). If VIV(7)=2., the state varia- 
bles Vj., y, r, AZ, geodetic latitude 0 , and 
node w, must be input into VI V (1-6). 

If payload is required as a constraint 

VIVC8-10) should contain initial values 
of ideal velocity , back pressure loss, 
and gimbal loss, respectively* 


VRCUT 


FT/SEC Relative velocity at which the last thrust 
event of the first stage is to be cut off, 
if desired. Preset to zero. 


LBM/SEC Flowrate per engine for each thrust event. 

LBM Liftoff weight at TZER 0. 

LBM Maximum critical propellant in a stage from 

which performance reserves are taken. This 
value Is assumed to include the performance 
reserves and must be input if IPR is greater 
than zero. Preset to zero. 


WTBL 

(25) 

M/SEC 

Wind speeds corresponding to altitude in 
ALTTBL and aximuth angles in AZWTBL. 

WTJET 

(15) 

LBM 

Jettison weights per thrust event. The weight 
jettison occurs at the end of the thrust event, 
such as releasing the empty first stage. 

Preset to zero. 

WTLAND 


LB 

Landed weight of booster for flyback fueL 
calculation. 

XGP 

(15,2) 

M 

Longitudinal coordinate of engine thrust point 
location for each engine of each stage. See 
Figure 9-2. 

XLEN 

(2) 

M 

Aerodynamic reference length for each stage. 

XREF 


M 

Longitudinal coordinates of the aerodynamic 
reference point. See Figure 9-2. 

YGP 

(15,2) 

M 

Lateral coordinate of engine thrust point 
location for each engine of each stage. 
See Figure 9-2. 

ZGP 

05,2) 

M 

Vertical coordinate of engine thrust point 
locations for each engine of each stage. See 


Figure 9-2. 
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NAME 


DIMENSION 


UNITS 


EXPLANATION 


ZREF 


M Vertical coordinate of the aerodynamic moment 

reference point. See Figure 9-2. 
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Table 9-1. CONSTRAINT CODES 


The codes contained in this table are input into KCDPHT and KCDRES to 
designate the pwoff nnd the Interred \ ate and terminal constraints desired 
for the trajectory. The appropriate values desired for these constraints must 
then he input into PS1RF.0 and PSTRST. 


CODE NUMBER 

UNITS 

CONSTRAINT OR PAYOFF 




1 

KG 

Payload 

2 

M/SEC 

Inertial velocity 

3 

DEG. 

Inertial flight path angle 

4 

M ? 7 

Radius 

5 

6 

Mf/SEC 4 

tfVSEC 

Energy 

Angular momentum 

7 

DEG. 

Inertial longitude 

8 

DEG. 

Inertial heading angle (+ East from South) 

9 

DEG. 

Colatitude 

10 

DEG. 

Inclination 

11 

DEG. 

Line of modes 

12 

M 

Semi-latu s rectum 

13 


Eccentricity 

14 

SEC 

Bum Time 

15 

LB/FT 

Ma;:itnum dynamic pressure 

16 

DEG. 

True anomaly 

17 

DEG. 

Argument of perigee 

18 

DEG. 

Not used 

19 


Not used 

20 

NM 

Flyback range 


Table 9-2. OPTIMIZATION CODING 

The KDB locations listed below will provide optimization of the corresponding 
parameter when the location contains a 1. For no optimization of given para- 
meter, the parameter is either held at a constant value or is determined as a 
consequency of other parameters. All KBD locations contain zero unless input 
otherwise. 


KDB LOCATION 

UNITS 

CONTROL PARAMETER 

1-10 

SEC 

TAUT (1- ID) ( burn 1 !»..•> 

1 ) 

KG 

W01 (liftoff weight) 

12 

DEG 

AZ (launch azimuth) 

13 

DEG 

CPTBL(2) Booster pitch attitude at time 



TTBL(2) 

14 

DEG 

CPTBL(3) Booster pit ;h attitude at time 



TTBL(3) 

15 

DEG 

CPTBL(4) Booster pitch attitude at time 



TTBL(4) 

16 

DEG 

CPTBL(5) Booster pitch attitude at time 



t TTBL<5) 
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17 

DEG 

CYTBL(2) Booster yaw attitude at time 
TTBL(12) 

18 

DEG 

CYTBL(3) Booster yaw attitude at time 
TTBL(n) 

19 

DEG 

CP0TBL(l)Orbiter pitch attitude at time 
T0BL(1) 

20 

DEG 

CP0TBL(2)Orbiter pitch attitude at time 
T0BL(2) 

21 

DEG 

CP0TBL(3)Orbiter pitch attitude at time 
T0BL(3) 

22 

DEG 

CY0TBL(l)Orbiter yaw attitude at time 
T0BLU1) 

23 

DEG 

CYT0BL(2)Orblter yaw attitude at time 
T0RBU2) 


Table 9-3. THRUST EVENT CUTOFF OPTIONS 

The table below lists the triggers that are turned on for a particular 
thrust event according to the value of MSWCH(ITHR' V , where ITHR Is the number 
of the thrust event. The thrust event will be terminated when the trajectory 
reaches the first trigger value selected for that event, except for the g-limit 
trigger when used with continuous throttling, in which case the trigger 
(if it Is the first trigger reached during the thrust event) will cause only a 
print of the state at the time when the g-limit was reached. The thrust 
event will continue then until the next trigger is reached, with the thrust 
adjusted so that the g-limit is held for the remainder of the thrust event. 

The absolute value of MSWCH is used to denote the triggers to be turned on, 
while the algebraic sign tells the program whether to use continuous (+) or 
discrete (-) throttling at the g-limit. 


ABS (hSWCH(ITHR) ) 

TRIGGER USED 

0 

TIME 

1 

TIME, G-LIMIT 

2 

TIME, VSUBR 

3 

TIME, G-LIMIT, VSUBR 

4 

FUEL 

5 

FUEL, G-LIMIT 

6 

FUEL, VSUBR 

7 

FUEL, G-LIMIT, VSUBR 
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Table 9-4. USE OF LPRINT OPTION 


BASIS FOR 
PARAMETER 
UPOATE 

LPRINT = 2 

LPRINT - 1 

LPRINT = 0 

Constraint 

Restoration 

Step 

Trajectory Block Output 

The trajectory block Is 
printed for every iteration*. 

Trajectory Block Output 

Only first and final tra- 
jectory blocks are printed. 
The weight summary table is 
provided at the end of the 
restoration step printout. 

Trajectory Block Output 
Same as LPRINT • 1 

Optimization Procedure Output 
Output formats 1 thru 6 

Optimization Procedure Output 
Same as LPRINT = 2 

Optimization Procedure Output 
Output formats 1.5, and 6 

PAYOFF 

OPTIMIZATION 

STEP 

Trajectory Block Output 
Same as above 

Trajectory Block Output 

Only first and final tra- 
jectory blocks are printed. 
The weight summary table Is 
printed at the end of the 
optimization step. 

Trajectory Block Output 
Same as LPRINT « 1 

Optimization Procedure Output 
Output formats 1 thru 8 

Optimization Procedure Output 
Same as LPRINT = 2 

Optimization Procedure Output 
Output formats 1 . 5, 6, 7, & 8 


4 An iteration is defined as a pemanant parameter update. 







9.3 DETAILED OUTPUT DESCRIPTION 

The RAGMjflP output is divided into two groups; 1) the trajectory output, 
which displays the running account of the restoration and optimization pro- 
cedure, and 2) postprocessor output in the form of publishable output tables 
and plots* This subsection describes the output of the trajectory and the 
input requirements and the output of the postprocessor modules. 


9.3.) Trajectory Output 

The RAGM0P trajectory output can be classified into the following three 
general areas : 


• Namelist Input , All Namelist data is printed describing the current 
status of each variable in the Namelist input array* Included with 
the Namelist, immediately preceding the first Namelist variable, 
the locations(in meters) of the equivalent engines used in the moment 
balance scheme for each stage are printed for reference* 


• Trajectory Block Output . Variation of the equations of motion as well 
as other pertinent data is printed by a block printout based on 
print interval, discrete events, or discontinuities in the equations 
of motion. A page of output consists of a header row printed at the 
top of % e page (containing the case number, the HEAD input card 
holleritn field, and the page number) and four block printouts (see 
Figure 9*3). If the block format is describing an event, the block 
will have displayed to its left one of the following captions; 

LIFT0FF - initiation of the trajectory from the launch site. 

IGNITI0N - jump start. 

THRUST^EVENT - specifies the beginning of a thrust event. 

END. TILT - specifies beginning of angle-of-attack control history 
used only if FACT>0). 

10.KMS - 10 kilometers altitude has been reached. 

14 a KMS - 14 kilometers altitude has been reached. 

(^MAXIMUM - maximum value of dynamic pressure* 

MACHINE - the value MACH « 1 has been reached. 

BEGIN A GLIMIT - specifies the beginning of either a discrete or continuous 
throttling event. 

INJECTION - intermediate and terminal trajectory constraints are imposed 
at these time points. 
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The block printout consists o£ eleven lines of output having the following 
labeled mnemonics: 


Line 1 


TIME - time from liftoff (sec) 

RANGE - relative range from the launch site to vehicle subpoint (n mi) 

RNGAN - inertial range angle between the launch site and the vehicle 
radius vector (deg) 

THRST - total vehicle thrust (lbs) 

XMLB - total weight (lbs) 

LACC - longitudinal acceleration (g's) 


Line 2 


Z13-X - Apollo 13 Z position coordinate, RAGM0P X coordinate 
X13-Y - Apollo 13 X position coordinate, RAGM0P Y coordinate 
Y13-Z - Apollo 13 Y position coordinate, RAGM0P Z coordinate 
ZU13X - Apollo 13 Z velocity coordinate, RAGM0P X coordinate 
XD13Y - Apollo 13 X velocity coordinate, RAGM0P Y coordinate 
YD13Z - Apollo 13 Y velocity coordinate, RAGM0P Z coordinate 


(m). 

(m). 

(o). 

(m/sec) . 
(m/sec) . 
(m/8ec) . 


Line 3 

R - radius from center of earth (m) 
VSUBI - inertial velocity (m/sec) 

GAMI - inertial flight path angle (deg) 
AZI - inertial azimuth (deg) 

LAT - geocentric latitude (deg) 

LONG - relative longitude (deg) 

Line 4 

ALT - altitude above ellipsoid (m) . 
VSUBR - relative velocity (m/sec). 

GAMR - relative flight path angle (deg) 
AZR - relative azimuth (deg) 

N0DS - inertial descending node (deg) 
INCL - inclination (deg) 


9-23 


Line 5 


GDLAT - geodetic latitude (deg) 

LUMP - latitude of the instantaneous impact point (deg) 
LNGMP - longitude of the instantaneous impact point (deg) 
CHIP - pitch attitude (deg) 

CHIY - yaw attitude (deg) 

CHIR - roll attitude (deg) 

Line 6 

MACH - Mach number (unitless) 

ALPHA - angle of attack (deg) 

BETA - sideslip angle (deg) 

DELPC - pitch global command (deg) 

DELYC - yaw girbal command (deg) 

DELRC - roll glmbal command (set=0) (deg) 

Line 7 

2 

Q - dynamic pressure (lbs/ft ) 

QALPH - product of dynamic pressure and 
QBETA - product of dynamic pressure and 
PIT M - total pitch moment (newton-m) 

YAW M - total yaw moment (newton-m) 

R0L M - total roll moment (newton-m) 

Line 8 

FAA - aerodynamic axial force (lbs) 

FAS - aerodynamic side force (lbs) 

FAN - aerodynamic normal force (1M) 

WD0T - propellant flowrate of liquid engines (lbs/sec) 
T0RB - thrust of the orbiter engines (lbs) 

DW0RB - weight loss of the orbiter engines (lbs) 


2 

angle of attack (lbs deg/ft ) 

2 

sideslip angle (lbs deg/ ft ) 
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QCN - product-of dynamic pressure and the aerodynamic normal force coefficient 
(Ibe/fO 

WIND - velocity of the input wind (m/sec) 

AZW • direction to which the wind is blowing (deg) 

ISP - liquid engine specific Impulse (see) 

TSRM - thrust of the SRM engines (lbs) 

DWSRM - weight loss of the SRM engines (lbs) 


Line 10 


CH.VL - characteristic velocity (m/sec) 
TN„l - turning loss (m/sec) 

GV„L - gravitational loss (m/sec) 

DR„L - drag loss (m/sec) 

BKP,L - back pressure loss (m/sec) 
G1M.L - gimbal loss (m/sec) 


(See Appendix F) 


Line 11 

ID^VL - ideal velocity (see Appendix F) 

GIMAN - pitch gimbal angle of two engine equivalent for zero aerodynamic 
, or SRM moments 

I 

i 

* 

0 

I 

I 

*1 

f 

t 
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Figure 9-3. TRAJECTORY BLOCK PRINTOUT 


jg REPRODUCIBILITY of THE ORIGINAL page is poor. || 



After completion of the trajectory simulation, a second block of data is 
output which presents a weight summary of the subsystems of the vehicle. 

These weights are used in the definition of the payload weight. This block is 
shown in Figure 9-4. The variables are defined as follows: 

WP0RB - orblter propellant burned (lb) 

CHVEL - characteristic velocity (ft/sec) 

DELVP - delta velocity for performance reserves (ft/sec) 

FPR - flight performance reserves (lb) 

WRESID - computed residual weight (lb) 

CRSDWT - constant residual weight (lb) 

WTANK - computed tank weight (lb) 

CTNKWT - constant tank weight (lb) 

WDR0P - drop weight (lb) 

WP0MS - OMS propellant weigh, (lb) 

PAYL0D - payload (lb) (the value will be negative when payload is the 

payoff since maximisation is performed by minimizing a negative 
number.) 
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Figure 9-4. PAYLOAD HEIGHT BLOCK 


A trajectory constraint summary table Is then printed which presents the 
payoff, the error in the desired equality constraints to the present and 
previous trajectories, the amount of change requested from the previous run, 
the amount of change obtained, and the percentage of predictability (see 
Figure 9-5) . 
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Figure 9-5. TRAJECTORY CONSTRAINT SUMMARY TABLE 
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Optimization Procedure Output . Following the trajectory output are 
several variables, vectors, and matrices used in the optimization procedure. 

The input variable LPRINT however will limit the amount of printout obtained 
from this area. Table 9-4 lists the output obtained using the various LPRINT 
options. To separate the various outputs a number will be used which will 
agree with the nomenclative of Table 9-4. 

1) Influence Coefficients - This block of output provides the user with 
the various partial derivatives of the constraints with respect to 
the optimization parameters. The size of this output block is a 
function of the number of equality constraints (left to right on the 
page) and of the number of optimization parameters (top to bottom) . 

The header "F0RWARD DIFFERENCES" or "CENTRAL DIFFERENCES" relates 
to the method used for determining these numerical partial derivations 
(see Figure 9-6). 
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Figure 9-6. INFLUENCE COEFFICIENTS TABLE 


2) Weighting Matrix (HWIBT) - Identifies the weighting matrix used in 
the optimization procedure to normalize the various parameters 
(consult Section III) (See Figure 9-7). 
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Figure 9-7. WEIGHTING MATRIX OUTPUT 


3) Total I-SY-SY Matrix - Matrix products used in the optimization 
procedure (consult Section III) (see Figure 9-B) . 
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Figure 9-8. TOTAL I-SY-SY MATRIX OUTPUT 
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4) DRH0, WDS, GNU - The variables are defined in the following manner 
(see Figure 9-9) : 

DRH0 - equality constraint errors used in parameter update equation. 
WDS - product of total I-SY-SY matrix and DRH0 vector. 

GNU - matrix product used in optimization procedure. 

DftM0««D& ( 6NU 

• •6«*23>H4*un 

• “•» l010V«fU) 

• 2»4tjSSb»0) tisn;2<t2*«0] •)t»l«2^*0tl 

Figure 9-9. DRHO, WDS, GNU OPTIMIZATION PARAMETERS 

5) Cl, E2 - Values of the restoration and optimization control step 
variables E2 and El (see Figure 9-10) 



Figure 9-10. CONTROL STEP VARIABLES PRINTOUT 


6) Parameter update and convergence test block - The following variables 
are printed in this block (see Figure 9-11) . 

DPAR - optimization parameter update requirement (one per parameter) . 

0LD„PC0N - post convergence factor for parameters (one per parameter). 

NEW.PC0N - present convergence factor for parameter (one per 
parameter) . 

P. SCALE - convergence parameter normalization factor (one per 
parameter) . 

TAUT(I) - if time is being optimized, prints out updated time. 

W01 - if liftoff weight is being optimized, prints out updated weigjht 

AZ - if launch azimuth is being optimized, prints out updated launch 
azimuth. 

CPTBL(l) - if the booster pitch attitude is being optimized, prints 
out updated values. 

CYTBL(I) - if the booster yaw attitude is being optimized, prints 
out updated values. 

CP0TBL(1) - if the orblter pitch attitude is being optimized, prints 
out updated values. 

CY0TBL(1) - if the orbiter yaw attitude is being optimized, prints 
out updated values. 

BETC0N - convergence indicator (BETC0N*T - all convergence criteria 
have met in subroutine ANEWCH , BET C0N*F- nonconvergence) 
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Figure 9-11. PARAMETER UPDATE AND CONVERGENCE TEST BLOCK 


7) The QY (or El) search output - after the program has met the terminal 
conditions » the optimisation procedure starts by varying the value 

of the control step El (Section III refers to this parameter as QY). 
When performing this search several additional parameters are printed 
(see Figure 9-12) . 

T0TSLP - the slope of the performance index with respect to the 
control step variable El. 

NSRCH - indicates the number of points in the QY search. 
C0MP0SITE.PAY0FF INDEX - total derivative of the constrained payoff. 

QY - value of the control step used. 

PREDICTED PAY0FF INDEX - the estimated value of the payoff index using 

either the coefficients of the quadratic or 
cubic fits. 

QUADRATIC .FIT „F0R„QY - value of QY predicted by the quadratic 

approximation. 

CUBIC. FIT. F0R.QY - value of QY predicted by the cubic approximation. 

8) After the converged trajectory has been printed, the following output 
will occur: 

Orbital element summary table - the following injection varia -les are 

printed (See Figure 9-13) . 

TIME - time (sec). 

RANGE - relative range (m) . 

RANGE. ANGLE - relative range angle (deg). 

INERTIAL VEL0CITY - inertial velocity (m/sec) . 

RADIUS - radius (m) . 

FLIGHT. PATH .ANGLE - inertial flight-path angle (deg). 

1NCL - inclination (deg). 
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Figure 9-12. QY SEARCH OUTPUT 
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Figure 9-13. ORBITAL ELEMENT SUMMARY TABLE 
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DES.N0DE - descending node (deg). 

FLIGHT AZIMUTH - launch azimuth (deg) . 

GDLAT - geodetic latitude (deg) . 

GCLAi' - geocentric latitude (deg) . 

L0NG - longitude (deg). 

INERTIAL. AZIMUTH - inertial azimuth (deg) . 

2 

Cl - angular momentum magnitude (m / sec). 

WEIGHT - mass (lb). 

2 2 

C3 - vis viva integral (twice the energy) (m /sec ) . 

ECNTRICY - orbital eccentricity (unitless) . 

AP0GEE. RADIUS - radius of apogee (m) . 

PERIGEE. RADIUS - radius of perigee (m). 

AP0GEE. ALTITUDE - apogee altitude above ellipisoid (n mi) . 

PERIGEE. ALTITUDE - perigee altitude above ellipisoid (n mi) . 

After the orbital element summary table the parameter set that yields the 
converged trajectory is printed (see Figure 9-14). 

Parameter summary table — The values of all parameters, whether 

optimized or fixed, are printed for reference. 









TftUl 








r aflAf)pnapft r>2 

. 53000000*02 

.640»iS40*02 

•13191410*02 

•JiQISttO’QJ 

tOBflfflDQfltt r 000110000 

•00000000 


,00000000 

.rooocccc 

•ococcooc 

•cocoocco 

•ooooocoo 

.00000000 .00000000 


noi 

.5SCSCC00 *0 1 







II 
















- ? • r 1 7L04*2 1 

,*9JC«B19*CI 

.14BLS6C J«02 

•OLOOOOOO 

•CCOOOCCO 



: t\9t 

-.41780 Ltt-GS 

.oooocccc 

.cooorcco 

•CCOOOCCO 

•coooooco 



rwlfl 


• 9882 ?323*C2 

.oroccoor 

•OLOOOOOO 

•ooooooco 




•OOOOOOOO 

•OCOCOOQG 

•ooooccor 

•0000 00 OG 

•00000000 




Figure 9-14. CHI PARAMETER VALUES 


The last printout of the converged trajectory is the comment, 


"TELL EM RAGMOP DID IT". 
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9.3.2 Postprocessor Option 

The user has the option using the variable NTABLE to require the execu- 
tion of the postprocessor nodules B0PPTL and B0PTBL. The following paragraphs 
describe both the input requirements and the output for use of these modules. 

9.3. 2.1 B0PTBL Input /Output . The selection of the first or third option 
under the input integer variable NTABLE indicates the use of the output report 
table and subroutine, B0PTBL. A second Namelist input package, $INFUT2, is 
required after the trajectory input, $INPUT. The definitions for this 
Namelist package are described in Table 9-5. 

Table 9-5. B0PTBL INPUT VARIABLES 


VARIABLE 

DIMENSION 

DEFINITION 

PRESET VALUE 

DATE 

2* 

Hollerith Field for calendar data 

BLANK 

NCASE 

1 

Case number desired on output 

0 

OFFICE 

2* 

Hollerith Field for Office 
Identification 

BLANK 

SRID 

60* 

Hollerith Field for optional thrust 
event notation. One notation 
requires 2 fields. (Replaces the 
header "THRUST EVENT" with data found 
in SRID) 

BLANK 

TITLE 

8* 

Hollerith Field for desired table 

BLANK 


*A six letter Hollerith field makes up one dimension length . 


Output for subroutine B0PTBL consists of fifteen tables suitable for 
publication. The first seven tables give the output variables in the MRS sys- 
tem, the next seven tables give the same variables in the English system, and 
the last table is a table of contents and a list of definitions for the varia- 
bles in the various tables. Presented in Figure 9-15 are examples of the out- 
put tables one through fourteen. Table fifteen is shown first to describe the 
location and definition of the variables. Also shown is a data setup for the 
sample problem. 
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Figure 9-15. BOPTBL OUTPUT 
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Figure 9-15. BOPTBL OUTPUT (Continued) 
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Figure 9-15. BOPTBL OUTPUT (Concluded) 
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9. 3. 2. 2. BOPPLT Input /Out put . The selection of the second or third ootion un- 
der integer variable NTABLt indicates the use of the output CALCOMP plotting 
subroutine, B0PPLT. If the second option is used, the B0PPLT subroutine is 
used and required input for this subroutine directly follows the trajectory 
input; however, if the third option is used, the B0PPLT input follows the 
B0PTBL input. The input for the B0PPLT subroutine consists of a series of 
formatted input cards. These formatted input cards wi.il be described as a 
series of card sets. 

Card Set 1 . (13A6/13A6/A6) The first card set consists of three cards making 

up 27 Hollerith fields. (The first two cards consist of 13 hollerith fields 
each and the third card has one hollerith field, all start in column 1). This 
card set is used for one-way communication between the user and the operator 
of the CALCOMP plotter. (An example of this communication would be to tell 
the operator where to send the plot, what pen to use for the plotter, and 
what color ink to use.) 

Card Set 2 . (4A6,2E10.5,4A6,E10.5,I2) This card set will be referred to as 

the plot card since the information on th • card is used to determine what will 
be plotted. 

The input card is made up of a combination of alphanumeric and numeric 
words which are used in the description of the required plot. Alphanumeric 
words on the plot card are used in the label .ig of the axis on the CALCOMP 
plotter. Mnemonics corresponding to this card set are INAME(l-4) , XFAC, 

XMAX, INAME(5-8), YFAC, N and are defined as follows: 



INAME(l) 

- the alphanumeric name of the X variable required for plotting 
(see Table 9*6 for option^) 


INAMK(2) 

- the alphanumeric name of the units of the X variable required 
for plotting (see Table 9*7 for options) 

j 

INAME ( J-4 ) 

- additional alphanumeric words required by user on X axis label 
(optional) 


XFAC 

input multiplication factor (optional, if less than or equal to 
zero, XFAC is set to one) 

* 

X MAX 

* The maximum value of the X component required for plotting 
(optional, set to l.E + 20 if input less than or equal to 
zero) 

* 

\ 

I 

INAME (5) 

- the alphanumeric name of the Y variable required for plotting 
(s^e Table 9*6 for options) 
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INAME(6) 

lNAME(7-8) 

N 


the alphanumeric name of the units of the Y variable required 
for plotting (see Table 9-7 for options) 

additional alphanumeric words required by user on Y axis label, 
(optional) 

the number of label cards to follow the plot used. 


Table 9-6. ALPHANUMERIC PARAMETER WORDS FOR B0PPLT INPUT 


ALPHA 

CHIR 

DWSRM 

GIM.L 

L0NG 

QCN 

T0RB 

XD0T 

ALT 

CHIY 

FAA 

GV a L 

LTIMP 

RADIUS 

TSRM 

Y 

AZI 

CH,VL 

FAN 

ID a VL 

MACH 

RANGE 

VI 

YAW.M 

AZR 

DELPC 

FAS 

INCL 

N0DE 

RNGAN 

VR 

YD0T 

AZW 

DELRC 

GAMMAI 

ISP 

PIT.M 

R0L A M 

VWIND 

Z 

BETA 

DELYC 

GAMAR 

LACC 

Q 

THRUST 

WD0T 

ZD0T 

BKP a L 

DR A L 

GDLAT 

LAT 

QALPHA 

TIME 

WEIGHT 


CHIP 

DW0RB 

GIMAN 

LNGMP 

QBETA 

TN,L 

X 



See General Program Output , Section 9.3 for definitions of 
parameters. 


Table 9-7. ALPHANUMERIC UNIT WORDS FOR B0PPLT INPUT 


INPUT WORD 

DEFINITION 

INPUT WORD 

DEFINITION 


UNITLESS 


LB*DEG/FT**2 

BTU/FT 

BTU/FT**2 

METER 

METERS 

CAL/M* 

CAL/M**2 

M/SEC 

M/SEC 

DEGREE 

DEGREES 

NAUT M 

NAUTICAL MILES 

FEET 

FEET 

N*DEG 

N*DEG/M**2 

FT/SEC 

FT/SEC 

NEWTON 

NEWTONS 

G'S 

G’S 

N/M**2 

N/M**2 

KIL0GR 

KILOGRAM 

N-M/M* 

N-M/M**2 

KIL0ME 

KILOMETER 

POUNDS 

POUNDS 

LB/ FT* 
LB- FT/ 

LB/FT**2 

LB-FT/FT**2 

SECOND 

SECONDS 
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It is important that the user be careful in hia selection of the alpha- 
numeric words since the subroutine will only recognize the words of Tables 9-6 
and 9-7. If the user does misspell a word, the subroutine writes out the 

comment " IS NOT A RECOGNIZABLE PARAMETER WORD, CONSULT USER’S GUIDE/ 

THIS PLOT HAS BEEN ABORTED” if the word is a parameter, and IS NOT 

A RECOGNIZABLE UNIT WORD, CONSULT USER’S GUIDE/THIS PLOT HAS BEEN ABORTED” if 
the word is a unit. If either word is in error, the present plot is aborted 
and the next plot card is read. 

The plot card can also be used for two additional purposes; 1) a comment 
can be used for identification by placing a "C" in column ^.e with five spaces 
(see subsection 9.4.3). The rest of the card (other than the numeric portion) 
can be used for alphanumeric words. A comment card, when read is printed and 
ignored. 2) plot card can also be used to identify the end of the plot 
cards and the end of the plotting requirements. Placing the word "ENDPLT" 
starting in column one calls the final CALCOMP subroutine and terminates the 
plotting sequence (see subsection 9.4.3). 

Card Set 3 . (E°.5,2E1C.5,5A6) This card set, which wiJl be called the label 

card, is an optional set based on the value of N i^ lot card. If N=0, 

no label cards are required, and if N>0, the number label cards will be a 
function of the value of N. Mnemonic corresponding to the label card are the 
following: X, Y, YHT , L0CN where, 

X - the lateral location of the label referenced to the origin of the 
axes (input in inches) 

Y - the horizontal location of the label referenced to the origin of the 
axes (input in inches) 

YHT - the height of the letters In llu* label (Input In mitts <>1 1// o! an 

inch. For most labels a value of J is adequate). 

L0CN - Alphanumeric words used for the label (not to exceed 30 spaces). 

Output for the B0PPLT subroutine is shown in tigures 9-16 and 9-17. 
Figure 9-16 shows the paper output of the plot cards, whereas Figure 9-17 
shows the plotted information. 
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Figure 9-16. PAPER OUTPUT FOR B0PPLT 


WOLTOS 


BIN215 


43271 1 


RPCMOP 

CS0217 

072172 

075203 

Figure 9-17. 

B0PPLT PLOT OUTPUT 
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Figure 9-17. (Continued) 
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Figure 9-17. (Concluded) 
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9.4 CONSTRUCTION OP THE DATA DECK 


9.4.1 Control Cordi 

The card input sequence for a given case on the RAGMOP program is 
initiated by the computer system dependent control cards. The MSFC Univac 
1108 Exec "III system control cards are set up and arranged corresponding 
to the sample control data input of Figure 9-18. A brief description of the 
control card3 is presented below (for a better understanding of them or any 
other control cards the user should consu.’ c the appropriate Programmer's 
Procedures Manual . ) 



Figure 9-18. SAMPLE CONTROL SETUP FOR MSFC UNIVAC 1108 EXEC VIII SYSTOi 
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Card 1 . M v?RUN , / / T „RAGM0P , Seq # , NAMEB IN , T IME AGE/ CARDS 11 

This card is the run card which executes the control stream logic 
aud provides the computation laboratory with pertinent information 
concerning the name of the job, the job charge number, the name of 
the user, the user's bin number (for return), the maximum C.P.U. 
time of the run, the maximum number of pages, and the maximum 
number of output cards. 

Card 2 . ASG , T A RAGM0P , T , TAPE# " 

If the program is stored on magnetic tape, this card is used to 
identify the tape number and assign a name for the tape for 
execution purposes. 

Card 3 . ’ , (3PEWIND-RAGM0P ,, 

Assures a rewind of the program tapes to the first record. 
card 4 . n (?C0PIN.RAGM0P ,TPF$ 

Copies the program tc?e onto a temporary storage drum file. 

Card 5 . "CFREE.RAGM0F 

Frees the program tape from the control stream, rewinds the tape an-* 
turns off the tape unit- 

Card 6 . "@MAP , IS .M0P ,M0P ,, 

Processor call statement for the collector, the name tl M0P tl is given 
to the program file containing the binary elements of the program 
subroutines and the required external (library) subroutines. 

Card 7 . ” LIB, V SYS&*MSFC$." 

Control card to obtain MSFC library subroutine package. 

Card 8-J3 . Program overlay cards (required to get within 32k limits at MSFC) ; 
optional if the 32k case limit is not imposed.) 

Card 14 . "eXQTJ^P" 

Demands execution of the program file M0P. 

Trajectory input . Namelist input for trajectory simulation. 

M @F IN tT . Specifies finalization of the control flow. 


?.4.2 Namelist Input 

Program input (except Cor the B0PPLT input) uses the Namelist format. 
Namelist is a Forfar. V input routine which groups and stores the input into 
various sets of input data depending upon the alphanumeric name of the variables 


9-54 







f 

# 


and the dimension of the variables in the calling subroutine. Before using 

the Namelist type of input, the user should be aware of the following rules. 

1. No Namelist fnput (names and data) may start in column 1. 

2. The firs* control card of a Namelist group of data must contain a 
dollar sign ($) folic -^d by an identification name. This name is 
followed by at least one blank character. 

3. Namelist data must take one of the following forms 

a. variable name ® constant, where variable name may be an array 
elemeut name or a simple variable name. Subscripts must be integer 
constants . 

b. array name ® set of constants (separated by commas), where K* 
constant, may be included to represent K constants (K must be an 
unsigned integer). The following ire two examples of how TIM may 
be expressed. 

TIM - 52., 52., 52., 52., 

TIM - 4 * 52., 

c. subscript variable * set of constants (separated by commas). This 
causes elements to be loaded, starting with the element designated 
by the subscript, in consecutive order. 

4. All data elements must be separated by commas. 

5. Table 9-8 lists the types of constants which may be utilized with 
Namelist, the definition of the constants and their magnitude limits 
on the Univac 1108 digital computer. 


Table 9-8. NAMELIST CONSTANTS 

CONSTANT 

DEFINITION 

MAGNITUDE 

Integer 

1 to 11 digits written 
without a decimal point, 
(fixed point number) 

+ 10 10 

Real 

1 to 9 significant digits 
written with a decimal 
point (floating point number) 

10- 38 - 10 38 

Double 

Precision 

1 tD 18 significant aigits 
written with a decimal point 

CO 

o 

CO 

+ 

0 

1 

00 

o 

1 

o 


w 
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6. Fortran V assumes tha^ all variable names beginning with I, J, K» L, 

M or N are integers ana thi ** allows no decimal, (Unless the variable Is 
REAL in the input subroutine). Variable names beginning in other than 
these six lett^ must contain a decimal. 

7. Each card in ue Namelist which contains data must end with a comma. 

8. If an array is being input and the amount of data exceeds one card the 
data may be continued on the next card by either continuing the array 
or by redefining the variable. 

9. Real numbers may be written with or without the decimal exponent speci- 
fied. The following example shows three ways in which a real number may 
be written 

2230.5 
22305. E-l 
.22305E4 

A double precision constant must be followed by the letter D with a 
signed (+ is optional) one, two or three digit exponent. The following 
are acceptable double precision constants. 

1.0D0 (means 1.0) 

75.3D+1 
7 .5302 
75300. D-2 

10. The end of a Namelist requires a $END card either on the same card with 
the last data item or 33 the last card of the deck. 

11. The number of elements read into an array must not exceed the dimensions. 

12. Namelist input can include columns 2 through 80. 


The trajectory input for each case is initiated by a $INPUT card and 
terminated by a $END card. No particular order of input variables has to be 
maintained when setting up data cards, but it is suggested that alphabetizing 
the input or defining some order oi sequence should be used to minimize both 
checkout procedures and input errors. 


Multiple cases can be input by continued use of the $INPUT and $END cards. 
However if multiple cases are rec ired, tne variable LAS T n st be set equal to 
the number of cases in the initial data set to signify uhat subsequent cases 
will follow. Multiple case da 4 * a is made up of the variables des5 ed to be 
changed from the initial or preceding case. 



9.4.3 Input Cotugories 

Input to the RAGM0P program can be grouped into four general categories. 
These are 1) trajectory physical model, 2) program control, 3) trajectory 
optimization and control, and 4) postprocessor control. Tables 9-9 through 9-11 
define the variables used under each category. Definitions of these four 
categories are presented below. 

Trajectory Physical Model . This category relates to all variables concerned 
with either the earth model or the vehicle physical model. Subdivision of 
this category is used to classify various topics which raake up the physical 
model. These topics are the following. 

Earth Model Constants - constants associated with the earth model shape 
and gravity model (not required input, preset values are normally used). 

Earth Model Variables - variables associated with the earth model, such 
as launch condition (preset values can be used) . 

Vehicle Model Variables - variables associated with the definition of the 
vehicle model. These variables are broken into the areas of aerodynamic, 
vehicle geometry, flight performance reserves, flyback, payload, and propulsion 
system. 

Program Control . This category describes the input which i& critical for the 
setup and operation of the trajectory simulation. Included in this category 
dre the integration inputs, logic flags, and thrust event configuration re- 
quirements. 

Trajectory Optimization and Control . Thi9 category consists of all parameters 
associated with either the optimization or attitude control input requirements. 

Postprocessor Control .* The variable NT ABLE in the $INPUT Namelist is used to 
specify to the program that a postprocessing module will be initiated after 
completion of the trajectory simulation. 

* See eubparagra^ns 9. 3.2.1 and 9 .3.2.2 
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Table 9-9. RAGMC? NAMELIST INPUT SYMBOLS (ORDERED ACCORDING TO USE) - TRAJECTORY 
PHYSICAL MODEL 


INPUT SYMBOL 

DIMENSION 

BRIEF DEFINITION* 

PRESET VALUE 

UNITS | 

1 "TTH MODEL CONSTANTS 

CJ 

1 

Second ha**monTC in gravity potential equation 

. 1 62345E-03 

uni tless 

CMUE 

1 

Product of the universal gravitational constant and 
the Earth's mass 

T.986032E+14 

m 3 /sec 2 

DJ 

1 

Fourth narmonic in gravity potential equation 

.7875E-05 

unitless 

FLAT 

1 

Flattening coefficients of the Fischer ellipsoid 

1./298.3 

unitless 

G2ER0 

] 

Earth’s gravi tational acceleration constant 

9.80665 

m/sec 2 

H 


Third harmonic in gravity potential equation 

-5 . 75E-06 

unitless 

0MEGA 

’ 

Earth's angular rotation rate 

. 72921 158E-04 

rad/ sec 

RE 


Earth's equatorial radius 

6378165.0 

meters 

EAR T H MODEL VARIABLES 

ALAT 

1 

Geodetic latitude of the launch site (real variable) 

28.531855 

deg 

AL0NG 

1 

Longitude of the ’aunch site (measured positive west) 

80.5649528 

(j *>n 

ALTLS 

1 

Altitude of the launch site above model ellipsoid 

0. 

meters 

ALTTBL 

25 

Altitude table for wind and wind azimuth tables 

0. 

meters 

A2 

1 

Launch azimuth 

90. 

deg 

azwtbl 

25 

Wind azimuth table 

0. 

deg 

DTZ 

1 

Time from GKR to liftoff 

0. 

sec 

NWIND 

1 

Number of points in wind table 

0. 

uni tless 

WTBL 

25 

. 

Wind table 

. 

0. 

m/sec 

trajectory model 

VARIABLES 




Git" 

15 

Longi tudinal acceleration limit 

. 1 E20 

G's 

JTHR 

15 

- l # SRM are being used; - 0, no SRM are used 
(per thrust event) 



nswCH 

15 

Tnrust event control flag (see table 9.1) 

4 

unitless 

N0£V\T 

S 

Number of events per stage 

0 

uni tl ess 

M0MBAL 

T 

t 

■'oinent balance flag ( 0, no; % 1*2, 3, yes 

'see nan . 9.2); 

0 

um t less 


*Full description V3 given in Paragraph 9.3. 
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Table 9-9. (Continued) 


\?UJ SYMBOL 

DIMENSION 

BRIEF DEFINITION* 

PRESET VALUE 

UNITS 

T L I FT 

1 

Time indicating end of vertical rise 

8. 

sec 

TTIl.l 

1 

Time indicating end of tilt-over maneuver (optional) 

0. 

sec 

TZER0 

1 

Liftoff time 

0. 

sec 

VIV 

8 

Initial state vector for jump start 

0 . 

variable depending on value 
of ViV(T) 

I VEHICLE MODEL VARIABLES J 

| AERODYNAMICS 




CAALP 

(25.2) 

Coefficient of slope of CA w.r.t. angle-of-attack 
per stage 

0 . 

deg 

CA0 

(25.2) 

Coefficient of CA at zero angle-of-attack per stage 

0. 

unitless 

CLBETA 

(25,2) 

Coefficient of slope of CL w.r.t. sideslip angle per 
stage 

0. 

deg 

CMALP 

(?5.2) 

Coefficient of slope of CM w.r.t. angle of attack 
per stage 

0 . 

deg 

CD 

(25.2) 

Coefficient of CM at zero angle of attack per stage 

0 . 

unitU*ss 

CNALP 

(25.2) 

Coefficient of slope of "N (normal force) w.r.t. 
angle ot attack per stage 

0. 

deg 

CNBETA 

(25.2) 

Coefficient of slope of CN (yaw moment) w.r.t. 
sideslip angle per stage 

0. 

deg 

CNO 

(25,2) 

Coefficient of CN (normal force) at zero angle of 
attack per stage 

0 . 

unitless 

CYBETA 

(25,2/ 

Coefficient of slope of CY w.r.t. sideslip angle 
per stage 

0 . 

deg 

PNM 

(25,2) 

Mach number table for all aerodynamic tables & 
TALFP table pe» stage 

0. 

uni tless 

S 

15 

Aerodynamic reference area per thrust event 

0 . 

meters 

XLEN 

2 

Aerodynamic reference length per stage 

0. 

meters 

XREF 

2 

Longitudinal location of aerodynamic moment 
reference point per stage 

0. 

meters 

ZREF 

2 

Vertical location of aerodynamic moment reference 
point per stage 

0. 

meters 

BA .E AXIAL FORCE 

Ialtbas 

(25,2) 

Altitude table for base axial force table per stage 

0. 

meters 

Jbaxial 

(25,2) 

Base axial force table per stage 

0. 

lbs 


description is given in Paragraph 9.3. 




■Mmmd «M im 





9-60 


Table 9-9. (Continued) 


INPUT SYMBOL 

DIMENSION 

BRIET DEFINITION* 

PRESET VALUE 

UNITS 

GEOMETRY 




I* ONC 

. 1 

Reel pal of the nose cone radius 

0 . 

1/ft 

TBOWT 

(15.2) 

Delta weight history for center of gravity 
location per stage 

0 . 

lbs 

TXCG 

(15.2) 

Longitudinal component of the center of gravity 
per stage 

0 . 

meters 

TZCG 

(15.2) 

Vertical component of the center of gravity per 
stage 

0 . 

meters 

W 01 

1 

Liftoff weight 

0 . 

lbs 

WTJET 

15 

Jettison weights 

0 . 

lbs 

FLIGHT PERFORMANCE RESERVES (OPTIONAL METHOD) 

CWO0T 

15 

Cri tical flowrate 

0 . 

lbs/sec 

OELVG 

1 

iV for geometry reserves 

0 . 

m/sec 

OELVP 

1 

^V for performance reserves 

0 . 

m/sec 

WP* 

1 

Maximum cri tical propellant 

0 . 

lbs 

IPR 

1 

Specifies thrust event from which flight 
reserves are taken 

0 . 

uni t less 

FLYBACK 

3100 

1 

Booster lift to drag ratio 

0 . 

unitless 

GAFPCT 

] 

Percentage of go around fuel 

0 . 

decimal fraction 

E0PCT 

1 

Additional flyback fuel percentage for engine 
out (of loaded weight) 

0 . 

decimal traction 

HEDHND 

1 

Flyback headwind 

0 . 

ft/sec 

K.BACK 

1 

=1, flyback computation required, *0, no 
flyback 

0 . 

unitless 

SFC 

1 

Thrust specific fuel consumption of the flyback 
engines 

0. 

ft/lb-hr. 

VCROS 

1 

Flyback cruise velocity 

0 . 

ft/sec 

WTlAND 

1 

Landed weight of booster 

0 , 

lbs 


*Fu!l description is given in Paragraph 9.3. 


Table 9-9. (Concluded) 


| INPUT SYHBOL 

UIHENSION 

BRIEF DEFINITION 

PRESET VALUE 

[ UNITS 

| PAYLOAD COMPUTATION f 

C0RBWT 

1 

Constant orblter weight 

t 

lbs 

CRStIWT 

1 

Constant residual weight 

0. 

lbs. 

ctukwt 

1 

Constant tank weight 

0. 

lbs 

OV0HS 

1 

Velocity potential of OHS engine 

0. 

ft/sec 

FPRFAC 

1 

Scale factor for flight performance reserves 

0. 

unitless 

0MSISP 

1 

Specific Impulse of the QMS engine 

0. 

>ec 

SCALE 

1 

Scale factor for orblter tank 

0. 

unitless 

SCRSO 

1 

Residual scale factor 

0. 

uni tless 

VISPB 

1 

Specific Impulse for booster 

0. 

sec 

VISP0 

1 

Specific impulse for orblter 

0. 

sec 

1 PROPULSION SYSTEM I 

AE 

15 

Exit area per liquid engine per thrust event 

. 

0. 

■n 2 

CISP 

15 

Specific Impulse of liquid engines used for con- 
tinuous throttling per thrust event 

0 

sec 

flbs 

15 

Vacuctfn thrust per liquid engine per thrust event 

0. 

lbs 

FuELIQ 

15 

Liquid engines fuel per thrust event 

0. 

lbs. 

NPARBN 

1 

Parallel burn indicator. *1 , parallel burn liquid 

0. 

unitless 

SRHAE 

1 

SRH exit area per engine 

0. 

m 2 

SRDWTB 

15 

Total SRH weight overboard table 

0. 

lbs 

SRMfTB 

15 

SRH sea level thrust per engine table 

0. 

lbs 

$RH>RP 

1 

Total SRH fuel available 

0. 

lbs 

SRHTTB 

15 

Independent time table for SRDWTB & SRHFTB 

0. 

sec 

SRXCGf 

1 

Longitudinal SRH aiming point 

0. 

m 

SRXGP 

5 

Longitudinal SRH glmbal points 

0. 

HI 

SRYGP 

5 

Lateral SRH glmbal points 

0. 

m 

SR2CGF 

1 

Vertical SRH aiming point 

0. 

m 

SR2GP 

5 

Vertical SRH glutal points 

0. 

m 

THE 


Nuntoer of liquid engines per thrust event 
(TNE(1 .thrust event ))and number of SRH 

0. 

unitless 



engines oer thrust event (THE (2 .thrust event)) 



WO0T 

15 

(15,2) 

Propellant flowrate per liquid engine per thrust 
event 

0. 

lbs/sec 

AGP 

Longitudinal glmbal plan loca M on per Mquld 
engine per stage 

G. 

in 


YGP 

(15,2) 

Lateral glmbal plan location per liquid engine 
per stage 

0. 

m 

ZGP 

(15.2) 

1 

Lateral glmbal plan location per liquid engine 
per stage 

0. 

i m 


*Full description ie given in Paragraph 9.3. 
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Table 9-10. RAGMOP NAMELIST INPUT SYMBOLS (ORDERED ACCORDING TO USE) 
- PROGRAM CONTROL 


INPUT SYMBOL 

DIMENSION 

BRUT uEFlNlTUN* 

ratsn VALUE 

UNlls 

AVI 

' 1 

Integration error check verleble 

2 E-3 

uMt’ess 

i CASE 


c«te number 

1 

un1*'*ss 

1 EU 


Upotr bound on Integrity. error 

1 E-S 

unit less 

HE AO 

15 

Hollerith field for header output on trajectory 

BLANK 


MMN 

l 

Minimum step site for Integration 

25 

sec 


1 

Earth model f1eg, -0, oblate. *1 spherical 

0 

jnUless 

JU*» 

1 

Jur^ start flag 

1 

unitless 

KINO 

1 

Integration tyP* fla*. -1. var'able step-size 
Adems-MOult on; *2. Runge Kutta. »3, fined step 
Adams -Moti I ton ) 

1 

u« dess 

1 

KROER 


Order of differences In Integration 

i 

umtless 

LAST 


■0, If only one case is sum, *1 if more than one j 

case is ryn 

0 

j un'tiess 

IF*INT ! 

1 

Print indicator (see table 9-4} 

1 3 

; ur’t'.ess 

j t*Al 

1 

Hannwm number of Iterations 

I 0 

umtless 

1 NT ABLE 1 

1 

1 

r utcut table or plots required, *0, no tables or j 

plots. *1, tables only. *2 plots only, «3 tables | 

& plots j 

0 

unitless 

1 

i 

PRINT S 

15 

Desired print Interval per thrust event j 

10 

se* 

STEP 

^ 

Integration steo size per thrust event [ 

L i 

1 


Table 9-11. RAGMOP NAMELIST INPUT SYMBOLS (ORDERED ACCORDING TO USF.) - 
AND CONTROL 


INPUT SYMBOL 

DIMENSION 

BRIEF DEI INJTIQN* 

PRESET V»lcE 

— 

ALPHA 

1 

Attitude control flag, .optimize control, *1 , 

control based on zero aerodynamic normal forca, »2 , 
control based on ztro angle of attack, -3. control 
based on angle of attack history as a function of 
Mach number 

-1 

unifies, 

1 

CP0T8L 

20 

Pitch attitude history for orblter 

n 

‘ 

cptbl 

20 

pitch attitude history for booster 

iT 

deg 

rrBTBL 

20 

Yaw attitude history for orblter 

T 

i deg 

CYTBL 

20 

Yaw attitude history for booster 

3 

deg 

f N0 

10 

Tolerance on terminal end conditions 

0 

■ rariat'e 

KCDPHI 

10 

Terminal function codes (cede in KCDPHI (1) Is 
pay *f ) (see table 9 1 ) 

KCDPHI Mj-i 
KCDPMl (2/*2 
KCDPHI ( 3)* 3 
KCDPHI (4)»4 

^ umt’ess 

KCDRES 

6 

Intr mediate constraint codes (see table 9 1) 

n 

„r * 1 1 e . r 

kDB 

30 

Codes to identify optimization parameters 
(see table 9 2) 

0 

un’tiess 

kDT 

30 

Companion matn* to KDB for calculation of flight 
reference reserves 

o 

um t less 

KP 

* 

Order of **•.« deseed pitch polynomial per stage 

0 

ur • tless 

KV 

* 

Order of the desired yaw polynomial per stage 

0 

1 1 ess 

NVR5T 

1 

Intermediate constraint imoosed at termination of 
this thrust event number 

n 

uni t less 

tPSIREQ 

1 

t 

10 1 

| 

Required values of the terminal constraints at 
termination 

R5IRI H 1 ) - 
7H7A 11 
9 {»(( Ji2)*0 

P'jIRIUi D* 

m r n f(.,' 

variable 

jpMRM 

6 

Constraint valuet desired at restart point 

f, 

v#r I able 

JOALPM 

1 

Ma*lmum v«lu* of the product of dyne* 1 ' preilt/r* 
and angle of att*<k 

1 f 

ft deg/*' 

WTAM 

' 

M*»law*t vvlue rjf the produ* f of dynamic preisgra 
»nd Odes'lp angle 

! 1 *><) 

« dey/m J 

J5PPAW 

1 

P 1 1' h ati 1 tud* Inirpmeft Y f ,f me booster 

1 

-rey 

; V*AW 

1 

P I t r h attitude lorremenf fpr the '>r| ( iter 


deg 

I >Y*A0A 

1 

Y*w effffud* fnrreaaent for me bolster 

) 

deg 

>TAA» 1 

1 

1 

Yaw attitude increment for the orblter 

5 

1 '1*9 

Italfp | 

(».2> 

Ang> of attack control history per stage 

'l 

deg 

(taut j 

15 

Thrust event duration time per thryst event 

n 

! - 

i T0BL 1 

20 

Timetable for pitch A yaw attitude for orblter 

0 


TTBL ! 

20 

Timetable for pitch l yaw attitude for booster 

0 j 

| set 

VRCUT 

1 

Relative velocity cutoff 

0. 

j # T'S*r 


*Full description is given in Paragraph 9.3. 
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If NTABLE is equal to zero, no postprocessor is used and the data deck 
is set up corresponding to Figure 9-19; howt r, if NTABLE is greater than 
zero, three options are available which relate to NTABLE equal to one, two, 
or three. 

NTABLE equal to one signifies the use of the subroutine B0PTBL. This 
subroutine outputs the trajectory parameters into fifteen output tables 
suitable for publishing. The subroutine B0PTBL requires a separate group of 
NAMELIST input under the NAMELIST name $INPUT2. Therefore the input data set 
requires not only the $INPUT input but the $INPUT2 input (see Figure 9-20) . 

NTABLE equal to two signifies the use of subroutine B0PPLT. This sub- 
routine plots required variables using the CALCOMP plotter. Similar to 
B0PTBL, subroutine B0PPLT also requires a separate input package. However, 
formatted input is used rather than the NAMELIST form. When using the B0PPLT 
postprocessor the data package consists of the trajectory input plus the 
B0PPLT input requirements (see Figure 9-21) . 

NTABLE equal to three indicates that the subroutines B0PTBL and B0PPLT, 
respectively, are used for output tables and plots. In this option, the data 
is set up such that the B0PTBL input and B0PPLT input, respectively, follow 
the trajectory input (see Figure 9-22). 


'4INPUT 


TRAJECTORY INPUT (HTA8LE-0) 


CONTROL 


^RUN,//T 



Figure 9-19. SAMPLE DECK SETUP WITH NTABLE s 0 
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Figure 9-20. B0PTBL SAMPLE DECK SETUP 



I or HAX lARDS SHOULD B L INPUT 


Figure 9-21. B0PPLT SAMPLE DECK SETUP 
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figure 9-22. SAMPLE DECK SETUP USING B0PTBL AND B0PPLT 

9.5 SAMPLE PROBLEM 

The following pages contain excerpts from the computer printout of a 
typical sample problem. 

The vehicle used for this problem was the SRM shuttle with 049 HO orbiter. 
The mission was a 90-degree launch from Cape Kennedy to a 50 x 100-nautical 
mile orbit. Maximum payload for fixed liftoff weight was desired. Inter- 
mediate constraints are: 

qm = 3000 lb deg/ft 2 

8 llmlt “ 3 83 ( cont i nuo ' 8 throttling in the oibiter). 

I'arameterH aelectecl were; 

two values for a quadratic first stage x p polynomial, 

two x values for a linear second stage x polynomial with ri i «;coi.t I noons 
XI at staging, p 

and liftoff weight. 
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axor rag*o» 

FIRST BT»GE ratllTALENT EnGInF LOCATION c ARE 

I -.47955200*02, . 89746667 *0U. -•177? 9200 *02 1 T N IWut BOOT 


SECOND STAGE EBUIXAlENT ENGINE LOCATIONS ARE 

« .00000000 • .00000000 • .00000000 


* INPUT 
AE 


». »T 

along 

A„ T8 *S 


ALTLS 

A.TTl. 


ATL 

»? 

A 7WT9L 


I'.Xl ». 


• At 703 70GE*Ul, 
6 4 1 70 370OE* 01 . 
.000000006*00. 
W OOOOOOOOE* 00. 

• 286610606*0? 
.805649536*02 

. 0000000 06*00, 
. BS 34 40006.04 » 
.1 676400U6.G5, 
'• 3 35 28000E ♦ 05 • 
.000000006*00, 
.ooogooooc.oq, 

.0 00000006* UO. 

.000000006*00, 
. ooooooooe* oo, 
.000000006*00* 

•-OOOOOOOOE* 00* 
•00 0000006*00, 
. 00000000 £* 00 , 
.OOOOOOOOE.OO 
. OOOOOOOOE* UQ, 
•000000001*00. 
.ooooooooe*oo. 
.ooooooooe*oo. 

.OOOOOOC -♦ oil, 

. OOOOOOOOE * 00 » 
.000000006*00 
W20000000E -0? 

• 900000006*0? 

. 90000000E* 0? • 
.OOOOOOOOE* 00. 
0 OOOOOOOOE* 00. 
.OOOOOOOOE* 00. 
I. OOOOOOOOE *00, 
. OOOOOOOOE* UO. 
L OOOOOOOOE *00 
. 444800006*04, 

. 62 275 1 UOf • 06 , 

• 2 624 55006*06, 

. 17 79 3UU0E* 05 . 
.OOOOOOOOE* 00. 
•OOOOOOOOE *00. 
•OOOOOOOOE* 00. 
.000000006*00. 
.000000006*00, 
0000000006*00. 
.000000006*00. 
.000000006*00. 


.N1 7037006*01. 
.000000006*00. 
.000000006*00* 
•OOOOOOOOE* 00. 


. 304800006*04, 
. 8834 006*04, 
.213360006*05, 
.457200006*05, 
.000000006*00. 
.000000006*00. 
• Of. '*000006*00, 
.000000006*00. 
.000000006*00. 
.000000006*00. 
.000000006*00. 
.OOOOOOOOE *00. 
.000000006*00 

•100000006*05, 
.OOOOOOOOE* 00, 
•OOOOOOOOE *00. 
.000000006*00, 
.000000006*00, 
.000000006*00, 


.900000006*02, 

•000000006*00, 

.000000006*00. 

.000000006*00, 

.000000006*00, 

.000000006*00. 

.1556 90006 *06. 
.627199006*06, 
.137447006*06. 
.000000006*00. 
.000000006*00. 
.000000006*00. 
.000000006*00. 
.000000006*00, 
.000000006*00. 
.000000006*00, 
.000000006*00, 
.000000006-00. 


> IN INPUT BOOT 


.41703 7006*31. 
.OOOOOOOOE *oO, 
. OwOOOQQDOE *00, 
.OOOOOOOOE *U0 


.6096000 ?4. 

.914400006*04 , 
• 25 9080006 *05. 
.100000006*06 , 
.000000006*00, 
.000000006*00, 
.000000006*00, 
.OOOOOOOOE *00. 
.(00000006*00, 
.OOOOOOOOE *00. 
.OOOOOOOOE *00. 
.00000000F*OU. 


. lOOOOUQOE *0 7, 
.OGwO 00006 *00. 
.000000006*00, 
.000000006*00, 
.000000006*00, 
.oooooooor*oo. 


• 9000000U6 *o2 . 
.00000000 6 * 00 , 
.000000006*00, 
.000000006*00, 
.000000006*00, 
.000000006*00. 

. 33351 ?005*U6. 
.622751006*06, 
.622730006*05, 
.000000006*00, 
.000000006*00, 
.000000006*00, 
.000000006*00, 
.000000006*00. 
.00000000 6*00. 
•000000006*00, 
.OOOOOOOOE *00. 
.000000006*00, 


Figure 9-23. NAM" LI ST INPUT . 


axis system 


axis system 


.41 7037006*01, 
.OOOOOOOOE *00, 
.000000006*00, 


.-76 *>000 006 *04. 
. 121920006*05. 
/P8 95 60006*05. 
.OOOOOOOOE *00, 
.000000006*00. 
•OOOOOOOOE *00, 
.000000006*00, 
.0 0000 0006*00, 
.000000006*00, 
.OOOOOOOOE *00, 
.000000006*0(1. 
■ ( OOCOOliOE *00, 


*000000006*00, 
.OOOOOOOOE *00, 
VOODOO 006*00. 
• OOOCOOOOE *00, 
.000000006*00, 
.000000006 *00, 


.00000 0006 *00. 
.000000006*00, 
.000000006*00, 
.OOOOOQOr^QO. 
.OOOOOOOOt *00, 
.000000006*00, 

.489301 006*06, 
.45816 7006*06, 
.*355 86000 6*05, 
.OOOOOOOOE *00, 
.000000006*00, 
.000000006*00, 
.000000006*00, 
.000000006*00, 
.000000006*00, 
.000000006*00, 
.OOOOOOOOE *00. 
.000000106*00. 
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8 LOO 
CAALP 


CAO 


CASE 

CisP 


CJ 

C.^ET* 


C8ALP 


.ooooonoo6*oo. 

•OOOOOOOOE.OO 
. OOOOOOO OE* 00. 

•ooooooooe.oo. 
. 000000006 * 00 . 
• ooooooooe.oo. 
• 000000006. 00. 
•OOOOOOOOE.OO* 
• 000000006. 00. 
•ooooooooe.oo. 
.000000006*00. 
•oooooouoe«oo. 

•OOOOOOOOE.OO. 

• OOOOOOOOE.OO • 
.000000006.00. 
.820G00i JE-01 • 

• It 6000006. UQ. 
• 23AG0000C* 00* 
.21 *000006.00. 

• 180000006*00 • 
•000000006*00. 
.ooooooooe.oo. 
.oogoooqoe*oo, 

• 00 OOOOUOE* OU. 
.OOOOOOOOE.OO. 

• OObOOOUOE * 00 . 
.OOOOOOOOE.OO, 
•OOOOOOOOE.OO. 

•ioqoooooe.oi 

•8506SCUDE.0*. 
•A S06S00UE.U 3. 
•OOOCOOOOE.UU* 
•OOOOOOOOE.OO. 

• I6 23ASOUE -0? 
-r8 80000GliE-u2. 
-.Si OOOOUOE -0? • 
-.-S*000000E-o2. 
-.28 OOOOUOE -02. 
-. 1 *0000006- U% 

.00 OOOOUOE. 00# 
•OOOOOOOOE.OO. 
•OOOOOOOOE.OO. 
•QOOUOOOOE.UU. 
•OOUOOOUOE.OO. 
.OOOOOOOOE.OO. 

• OOOUOOUOE • 00 • 
•OOOOOOOOE.OO. 
-0 12000000E-01 . 
-. 135000006-01. 
-• 90000000E-0? . 
• SOCOOOOOE- 03. 
•200U0000E -01 • 
•OOOOOOOOE.OO. 
•OOOOOOOOE.OO. 
•OOOOOOOOE.OO. 
•OOOOOOOOE.OO. 


•OOOOOOOOE.OO 

•OOOOOOOOF *00. 
.000000006*00. 
.OOOOOOOOE.OO, 
.000000006*00. 
.000000006*00. 
.000000006*00. 
•000000006*00* 
.000000006*00. 
•OOOQOOOOE *00. 
.000000006*00. 
•000000006*00. 
.000000006*00. 
.000000006*00 
•82800C0CC-01 • 
.176000006*00, 
• 226000006*00. 
.201500006*00. 
•180000006*00, 
.OOOOOOOOF *00. 
. 000000006* 00 • 
.000000006*00. 
.000000006*00. 
.000000006*00. 
•000000006*00. 
•000000006*00. 
.000000006*00 

. *50650066*03. 
.*50650006*03. 
.003000006*00. 
•000000006*00. 

- .*80000006-02, 
-.580000006-02. 
-.500000006-02. 
-.230000006-02. 
- .180000006-02. 
.000000006*00. 
.000000006*00. 
.000000066*00, 
•000000006*00. 
•000000006*00. 
.000000006*00. 
•000000006*00. 
.000000006*00 
-. 12000000E -0 1 . 
- .2*0000006-01. 
-*660000006-02. 
• 1UUOOOOOE-02. 
.200 000 006 -01 • 
.000000006*00. 
.000000006*00. 
•OOOOOOOOF *00, 
.000000006*00. 


.000000006*00. 

.000000006*00. 

.000000006*00* 

.000000006*00* 

.000000006*00. 

•000000006*00. 

.000000006*00* 

.000000006*00* 

.000000006*00* 

.000000006*00* 

.000000006*00* 

.OOOOOOOOE*UO* 

.850000006-01* 
.250000006*00, 
.22*500006*00. 
.183500006*00* 
.000000006 *00* 
.000000006*00* 
.000000006*00. 
.000000006*00* 
.000000006*00* 
. 00000000 E *00. 
.000000006*00. 
.000000006*00* 


.*50650006*03 , 
.OOOOOOOOF *00. 
.000000006*00* 
.000000006*00 

- .*80000006-02* 
-.65000000r-u7* 
-.*50000006-02* 
-.200000006-02* 
.000000006*00* 
.000000006*00* 
•000000006*00. 
.000000006*00. 
.000000006*00. 
•000000006*00. 
.000000006*00. 
.000000006*00. 

-.12*000006-01. 
-.160000006-01. 
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Figure 9-23. (Concluded) 
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Figure 9-27. (Continued) 
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Figure 9-27. (Continued) 
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Figure S-27. (Continued) 
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Figure 9-27. (Continued) 
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Figure 9-27. (Concluded) 
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Appendix A 

COORDINATE TRANSFORMATIONS 


% 

\ 

* 

i 


A vector or tensor in any orthogonal coordinate system can be trans- 
formed into any other orthogonal coordinate system by premultiplication of 
the vector or tensor with the proper transformation matrix. Any coordinate 
system can be obtained from any other by at most three successive rotations 
about the coordinate axes. For any rotation about an x (or 1) axis, the 
transformation matrix is given by 

"1 


A , 
xx 


10 0 
0 cos 9 sin0 

0 -sin0 cos 9 


Aj, where x'y*z’ are 


the new axes. A vector in xyz is transformed into x'y'z* by premultiplication 
with A xx” e * 8 ‘ 

*’ ■ A xx’ i 


or 



1 0 

0 cos0 

0 -sin0 


0 

sinGl 

COS0 


where x' = x’e^, + y'e 2 » + z'e^' and x * xe^ + ye 2 + ze^ • 

Similarly, for any rotation through the angle $ about a y (or 2) axis the 
transformation matrix is: 


A x'x* 


COS<J> 

0 

-sin<f 

0 

1 

0 

sin<J> 

0 

COS if 


where x’y'z 1 is rotated about y* to form x*y*z*. Also, for a rotation of y 
about the z* axis we would have: 


* 

t 


A-l 


1 

* 

t 


« 





x*X 


cosy 

-siny 

0 


siny 0 

cosy 0 

0 1 


where the XYZ system is obtained by rotating x*y*z* about the z* the angl. y. 


Since any coordinate system can be obtained from any other by three suc- 
cessive rotations, the above matrices A^, A2* and A^ may be used together to 
form the transformation matrix from any one system to another. Note that for 
each. rotation, premultiplication of the proper matrix is required to transform 
the old coordinates into the new ones. Thus, if a coordinate system XYZ is 
obtained by rotating the x'y'z* system first about y an angle to form x*y*z*. 
then about x* an angle 0 to form xyz, and then about z an angle y to form 
XYZ, the transformation matrices are: 



cos4> 

0 


A , * = 

XX* 

0 

1 

0 


sin<|> 

0 

C03<P 


1 

0 

0 

A * = 
x*x 

0 

CO80 

slnd 


0 

-sin0 

CO80 




xX 


cosy 

siny 

o“ 

-siny 

cosy 

0 

_0 

0 



■ A„ 


and so the total transformation matrix is 


Vx 


A A . A » jl 
xX x*x x x* 


c 

1 -s 
0 


cosy 

siny 


siny 

cosy 

0 


10 0 
0 cos0 sin0 

|0 -sin0 cos0 


COS<t> 

0 

-sin<{> 

0 

1 

0 

sinij) 

0 

COS<t> 


A-2 





cosy 

sinycosS 

sinysine 


cos$ 

0 

-sin'Tj 

-siny 

cosy cos 9 

cosysinB 


0 

1 

0 

| 

0 

-sin0 

COS0 


sine)) 

0 

cos$ 1 


A 


x'X 


(co8YCos<H^inYsin9sin<t>) 

(-siuYCos^cosYsin0sln(|») 

cos9sin$ 


sinYcos© 
cosy cos0 
-sinO 


(-cosYsin^4-sinYsin9cos<ft] 

(sinYsin<|>+cosYSiP n ~os^ 

COS0COS<J> 


J 


Note that the order of rotations is important, buc .^at, if the proper order is 

maintained, the total transformation matrix may be obtained from any orthogonal 

system to another in the manner demonstrated above. The t ‘ans formation of any 

vector or tensor in the x'y' 2 ' system above to the XYZ system is accomplished 

then by simple premultiplication with A , : 

x x 

= A x'X R x' 


A- 3 


Appendix B 

INTERPOLATION METHODS 


Rapid, consistent, and accurate interpolation of tabulated data is 
essential to the operation of RAGMOP. To this end, two very efficient inter- 
polation schemes are used in the program. The linear interpolation subroutine, 
AMULG, is used for the determination of base pressure force as a function of 
altitude, wind direction and speed as a function of altitude, and center-of- 
gravity location as a function of total vehicle delta weight. A cubic spline 
interpolation routine, SPLINE, which ensures continuous first derivatives of 
the data across the data points, is used to find the aerodynamic coefficients 
as functions of Mach number and to find the atmospheric properties as functions 
of altitude in subroutine PRA63. Both interpolation routines keep track of 
the lower data point from the interval used in the previous call to the routine 
for each particular table. In this manner the routine is not required to 
repeatedly start from the beginning of the table to find the current interval 
containing the independent variable. Also, both routines will extrapolate 
beyond the range of the independent variable data set. 

The linear interpolation routine, AMULG, searches for the interval contain- 
ing the current value of the independent variable. If, of course, x « x^ (a 
data point), then y ® y^. If the independent data point is less than x^, 

it is assumed th&t the last valid data point has been passed and that extrapo- 
lation beyond that point will now be required. (Note that this allows the 
nunber of data points to be less than the maximum number of points available 
in each table). If x < x^, then m is chosen as hf^, (where m denotes the 
nunber of the base point of the interval) and backward extrapolation is per- 
formed. Once the base point, x^, is located, the general linear interpolation 
is: 


x - x 


y, 


m-l 


m-l 


m 


- x 


m-l 


(y_ - 

m 


The cubic spline interpolation routine, SPLINE, first locates the 
independent variable vf thin the central span of a four point group, i.e.. 


B-l 



x < x < x . , . If x < x, , w»l. If all of the available points in the 
independent variable table are filled, extrapolation will be performed when 
x > x ffi (where m is the number of locations in the table). If < x i+i» 

implying (as in AMULG) that the last valid data point has been passed, the 
routine sets x^ + £ “ x f+i + and sets the y^ + 2 va l ues 1° all the dependent 
variable tables equal to the y i+ ^ va * ue8< Thus, if fewer than the total number 
of available points are used in the tables, the scheme assumes that the depen- 
dent variables remain unchanged after the last valid input point. Having 
established the interval containing x (x < x < x ), the dependent variables 

m ntrl 

are found as: 


y± = a t X 3 + b ± X 2 + c ± X + y m . 


(B-l) 


The coefficients a^, b^, and c^ are derived below. 


3 - ^m*~l y m 

8i Vl" X m 


_ y m - Vl 

s = 

m x - x , 
m m-1 


_ y mf2 y m4-l 

Jp = V2 ' Vl 


If hfI. s is chosen such that 
m 


s + s 

“ or, in other words. 


28 - s. 


If m ■ n-1, s is chosen such that 

p 

S +3 

s ■ — — » or * n ot * ier words 


s = 2s - s , 

p 0) 
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Equation (B-l) can be rewritten as 


Ay l • y ®U ~ y m ” ***1 + bAx l + cAx l 


where 


Ax, ■ x _, i - x 
1 rofl a 


This yields, dividing by Ax^: 
2 

s » aAx^ + bAx^ + c. 


The spline assumption is 


iz 

dx 


s + s 


m 


and 


iz 

dx 


s + s 


3aAx^ + 2bAx^ + c » 


m*-l 


Note that this leaves the first derivatives continuous as m changes. 


Solving for a, b, and c we find: 


a * 


b 


(— r -8 + 8 - 0 


Ax^ 
s + s 


c • 


m 


X - x_ 


With R 


Ax 


Ax, X_.i ~ X 

1 nrt-i m 


we have y • y +((aR+b)R+c)Ax 

III 
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(3aR + 2b) R + c. 





Since a, b, and c are functions of the data points only (and are there- 
fore constant for a given interval of the independent variable data set) , 
their values are stored and reused without recalculation if interpolation 
within the sane interval is required more than once. 


l 
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Appendix C 

SOLUTION OF MOMENT BALANCE EQUATIONS 


From Paragraph 3. 2, 5. 2 we have the moment equations in the center-of- 
gravity/gimbal point coordinate system: 


(T tan 6 +T tan6 )dy + (T -T )dz = -Max 
y l y l y 2 y 2 y l y 2 

(C-l) 

(T tan 6 -T tan 6 )dz = -May 
y l P 1 y 2 P 2 

(C-2) 

(T tan5 +T tan 5 )dy «* -Maz 
y l p l y 2 p 2 

(C-3) 

6 « 6 * 6 
y l y 2 y 

(C-4) 


Since the total thrust of each engine is fixed 2) we can write: 


J 2 2 T ° 

T Z + T + T * Hr)" 
x i *1 *i 2 


(C-5) 


2 2 2 T 2 

and T + T + T - (-—) 
x 2 y 2 z 2 2 


(C-6) 


But, from the definitions of 6 and 6 : 

P y 


T » -T tanfi 

X 1 *1 Pi 


T m -T tan$. 

X 2 y 2 P 2 

T ”T tan6 «T tan6 

h *i y i H i 

T ■ T tan6 ■ T tan6 , 
z 2 y 2 y 2 y 2 y 


so that 


T/2 


y , j j— 

1 v 1 + tan 6 + tan 6 


(C-7) 


and 


T/2 


y , 2 j- 

1 + tan 6 + tan 6 


(C-8) 
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s 


Ue now assume that 6,6, and 6 are small enough so that 

Pi p 2 y 

2 

tan o << 1 

pi 

2 

tan 6 « 1 


and tan 6 « 1 

y 


(C-9) 

(c 10) 

(C-ll) 


This approximation yields, from equation (C-7) and (C-8), 

T 

T 

T * - - -- 

n 2 

T 

and T - 

y 2 2 . 

Using (04), (C-12), and (C-13) in (C-l), (C-2) , and (C-3) we have 


(C-12) 

(C-13) 


(T^ tan6y)dy * -Max 


(C-14) 


P2 


— (tan 6 - tan 6^ )dz • -May 

« P 2 

T 

-r~ (tan 6 + tan6 )dy » -Max 

2 p, P„ 


(C-15) 


(C-16) 


which can be rewritten as: 


tan6 


_ Max 
y ” T T dy 


tan6 


- cand - - 


and 


tan5 + tan6 


p 2 V z 

2Maz 

P 2 “ ' V y * 


(C— 17) 
(C-18) 

(C-19) 


Adding (C-18) and (C-19) we obtain 


. . 1 . May Mazv 

tan 6 ~ ( v ■* + j — ) 

P 1 t t dz dy 


(C-20) 
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Subtracting (C-18) from (019) yields 


tan <5 

P 2 


1 , Maz May * 

T t 1 dy " dz J 


( 021 ) 


Since the total thrust is inherently enforced in the solution, the error 
amounts to an error in the moments produced by the engine to balance the aero- 
dynamic moments. (It is possible to obtain an approximate solution that pro- 
duces no net unbalanced moment but changes the total thrust slightly. The 
solution used in RAGMDP does not do this.) The philosophy in RAGMOF is that 
the purpose of the moment balance scheme — to determine the performance of a 
vehicle based on realistic thrust vectoring — is satisfied acceptably when 
the moment on the vehicle has been reduced to a very small level without chang- 
ing the overall thrust of the vehicle. 

As expected, for no aerodynamic moment 6 ® <5 - 6 = 0. The thrust 

P 1 p 2 y 

vector of each engine is parallel to the center of gravity/ gimbal point y 
axis when aerodynamic moments are zero. 


Note that the moment and distance components in the above solution are in 
the XYZ (center-of-gravity/gimbal point) coordinate system. The moments must be 
transformed into this system from the body axis X'Y'Z' system in which they are 
computed, and the resulting thrust components must be transformed from the XYZ 
system back to the X'Y'Z' system in order to determine the total forces acting 
on the vehicle. The XYZ system is obtained from the X'Y'Z' system by rotating 
about Z' the angle (-p), where 


P 


-tan' 1 X GP 


Y' 

GP 



(C— 22 ) 


Thus, the transformation matrix (see Appendix A) from y'Y'Z' to XYZ is 
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)i] 


~ cos (-p ) 

sin(-p) 

0 ~ 


i* 


cosp 

-sinp 

o" 


r - 

V 


* 

-sin(-p) 

cos (-p ) 

0 


J* 

' - 

4>llp 

cosp 

0 

i 

y 

bJ 


_0 

0 

1 


k 


_o 

0 

1 


k' 

J 


and the reverse transformation matrix is then the transpose of the above, or: 


M 







i' 


cosp 

sinp 

o“ 


i 

W 

as 

-sinp 

cosp 

0 

< 

j 

U’J 

■ 

_0 

0 

1_ 




The aerodynamic moment* in the XYZ system are, then: 


M AX “ M AX' C0Sp + M AV' 8lnp 
m ay - m ax’ sinp + M ay , cosp 

m az " m az' 


(C-23) 


We define the distances 


DX » X' -X' =X' -X' 

GP X eg GP 2 eg 


01 • V V v V Y 'gp 2 -*' 


eg 


02 • Z GP 2 ' Z ' eg'” <z ' GPj" 2 ' eg’ 


so that 


cos p 


-1 DX 
p * tan DY 
DY 


V DX 2 + DY 2 


sin p 


DX 


V DX 2 + DY 2 


(C-24) 


and 
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(C— 25) 


dx ■ DX cos p - DY sin p = 0 


^X 2 + DY 2 


dy = Da sin p + DY cos p = 
dz ■ DZ 

Using (C-21), (C-22), and (C-23) in (C-14), (C-17) and (C-18) we obtain 


y * - tan 


-1 


6p^=-tan 


-1 


M DX M 

M AY' DY m ax' 


t t d *I 

fi+ S 

V 



M * M ' M ' DX 
AZ ”AY AX JX 

DY DZ ~ DYDZ 


t t(i + 

DY 


6 * -tan 

p 2 


-1 


M * M ' M avI DX 
AZ - AY + AY* 


DY DZ 


DYDZ 


t t ,i+wc 

M w 


Note that the gimbal angles in the X'Y'Z' body avis coordinate syste. are 
given by 

. , , -1 DX 

6 ' * 6 - tan — = 6_ - P 

Pi P! DY 


X _ -1 DX 

6 P , “ tan m " 6 - * p 


'l 

p 2 


and 


S' -6 

y y 


The total thrust components in the X'Y'Z’ body axis coordinate system are 
then found by noting that 


T ' 

* T 1 tan6 

X 1 

*1 

T ' 

* T ' tan6 

X 2 

CM 

T ' 

- T ' tan6 

Z i 

yi 
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where TXX, TYY, and TZZ are the total thrust 
coordinate system. 
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Appendix D 

PAYLOAD CALCULATION 


RAGMOP includes a payload calculation so that, with proper input, the 
payload as well as the orbiter (last stage) cutoff weight can be optimized. 
This calculation is based on a two stage space shuttle type vehicle. Payload 
is determined from the orbiter cutoff weight by including the requirements 
for: 

(1) flight performance reserves (FPR) for both stages based on the ratio 
of initial to final weight and the vacuum 1 for each stage, 

(2) orbiter tank weight based on total second stage fuel (including 
reserves) using an input scale factor, and 

(3) orbital maneuvering system (OMS) propellant requirements based on 
the delta velocity required and the I of the OMS, 


The calculations in this option are used at all times to determine the 
final mass printed with each trajectory summary. These calculations are per- 
formed as follows: 

Lj rx U — IJ 

P 0 CO 

cons ORB ORB 

W„ 

AV„ = 6__ characteristic velocity* 

P FPR 

FPR = W (1 - e' AV p /g o I SP ) 

CU ORB ° 


W TANK * n * W P 


+ FPR) 


cons 


W * W 

ON ORBIT CO 


- FPR - W, 


ORB 


TANK 


“p 0MS • “on ORBIT (l ' e ' SVoHsVsP OMS 


W PAYLOAD “ W 0N ORBIT ~ W P QMS “ W CONST QRB 


see Appendix F 



where 


FPR * flight performance reserves for both stages 

2 

= sea level gravitational acceleration (9.80665 m/sec") 

1 = booster (first stage) vacuum specific impulse 

SP B 

1 • orbiter (second stage) vacuum specific impulse 

SP 0 

^cp - orbital maneuvering system vacuum specific impulse 
OMS 

AV , = delta velocity required of orbital maneuvering svstem 
OMS 

AVp = delta velocity required for flight performance reserves 

= booster cutoff weight 
B 

W * orbiter cutoff weight 

C0 ORB 

W = constant orbiter weight (dry orbiter shell, no fuel or OMS propellant, 

CONSTQjyj no tank weight) 


W Q = booster liftoff weight (total vehicle) 

B 

W * orbiter initial weight (at staging) 

ORB 

W = on orbit weight (orbiter cutoff weight minus fuel, FPR, and tank weight) 

ON ORBIT 

Includes payload, OMS propellant, and constant orbiter weight. 

VL » fuel consumed by orbiter 

Peons 

orbital maneuvering system propellant weight 
payload weight 
orbite; fuel tank weight 

decimal fraction used in FPR calculation (typically .01) 
scale factor used in tank weight calculation (typically .08155) 


p 

OMS 

W PAYL0AD = 

W = 

TANK 

6 * 

FPR 

n * 


Note that the payload calculation is performed at all times and that if 
orbiter cutoff weight is desired rather tha:. payload, the user should omit the 
input for: <S FpR (FPRFAC) , n (SCALE), AV QMS (DV0MS) , and W C0NST (C0RBWT) . 
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Appendix E 

FLYBACK FUEL CALCULATION 

A flyback fuel calculation is included in RAGMOP in order to more real- 
istically optimize the performance of reusable flyback space shuttle boosters. 
This calculation is based on the Breguet* range equation and is performed as 
follows: 


The flyback range required is determined by subroutine ASIMP using a 
trivariant lookup. Flyback range is tabulated in terms of range, velocity, 
radius, and flight path angle at staging. The flyback fuel required in order 
to return to the launch site is the amount of fuel required to cover the fly- 
back range, plus reserves. The jettisoned booster weight at staging includes 
the flyback re' - d rements. 

The Bregi range equation gives the vehicle range as a function of 
lift over drag at * constant cruise velocity with a constant thrust specific 
fuel consumption for jet engine powered vehicles: 

. V <% . “o (E-l) 


where 


range 


V ■ cruise velocity (ideally velocity at best C^/C^) 
; » thrust specific fuel consumption 


— ■ best ratio of lift coefficient to drag coefficient for the vehicle 
i) (ideally at cruise velocity) 

W Q - weight at beginning of flyback leg 

and W 1 = weight at end of flyback leg. 

* Dormaaah, D. 0., Shirley, S.S. and Conolly, T.F., Airplane Aerodynamics, Pp 
3X3-351. Pitman Publishing Company, New York, 1967. 
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The flyback fuel required, W^, is given by 

R * c t CD 


W = W, (e v 01 -1) 
P 1 


(E-2) 


In order to include the effects of engine out and go-around fuel requirements, 
the approximation is made in RAGMOP that each may be expressed as an input de- 
cimal fraction of the landed weight of the booster. The equation for flyback 
fuel used in the program is, then: R r 

V CL /CD 

",-"l U + Pct E0 + Pc tGA )( C -1) 

where Pct„ » decimal fraction of booster landed weight required for engine out 
EO flyback 

and Pct_ “ decimal fraction of booster landed weight required for go-around 
capability at landing site. 
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Appendix F 

VELOCITY LOSS EQUATIONS 

The velocity loss equations result from the calculation of the change in 
the inertial velocity, V , during the time interval t - t, . 

i 3 0 

It is obvious that 


w - w 


Vjdt 


(F-l) 


Since 


V 2 » V • V 
I I I 


(F-2) 


it follows that 


V * V 
I I 


(F-3) 


Therefore 


ft, 


4V I - W - 


b y • V 
I I 


dt 


(F-4) 


Three different thrust levels are now defined. The first. T , is the total 

v 

resultant thrust calculated as the square root of the sum of the squares of the 
thrust components. The second, T g * is the scalar sum of the thrusts of all en- 
gines. The third, T , is the scalar sum of the vacuum thrust of all engines. 

s 


Now, adding and subtracting each of these divided by mass to the right hand 
side of (A) leads to 


je 

1 


F-l 







AV t - 


r t b f T T T T 
_ _51 + _s . _s 

I m m m m 


T vac vac V 

S + §■ ^ 4- 

mm m 




dt 


(F-5) 


Since T is the largesc of the three, we choose it to define ideal velo- 
vac 

s 

city so that everything else will represent losses. This gives 


AV ! = 


u 


vac T s vac T - T V • V T 
s , s . v s , I I 

+ + 4- 


m 


m 


dt (F-6) 


In RAGMOP 


V - - 
I ra 


(F-7) 


where 


F 

X 


cosx p cosx, r sinx p -sinx y sinx p 


T + T - F 4X , 

xs xx AN 

F 

y 

s 

sinx p cosx y cosx p -sinx y cosx p 


T + T - F aa 
yy ys AA 

.V 


0 sinx coax 

y y 

1 

T + SIDE 

ZZ 


(F-8) 


Taking the dot product of (7) and V^/V^ leads to 
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V I • V I 


g + A ( (T + T ) COS a. + (T + T ) cosa 2 + T cosa ) 
s y ~ v v xs xx 1 yy ys z 22 ^ 


- - (F,„co8a. + F. .coscu - SIDEcosa-) 
m AN 1 AA l j 


(F-9) 


cosot^ = 


cosa 2 5 


cosa^ = 


cosx p 

-sinx. 


cosXySinXp 

COSXyCOBXp 

. 8inx y 
’-sinXySinXp ' 

-8xnx cosx 

y p 


(F-10) 


If there is a roll angle, x » to be considered also, the variables in (10) 
are used to re-define coso^, cosa^, and cosa^ as 

ton, ° -nx r | 

r cosa. coso 2 eofiOjl ^ Icosc^ cosa 2 coho 3 1 10 10 (F-ll) 

L J L * Lsinv 0 cosx 


sinx r 0 cosx r 


Rewriting (6) and using either (10) or (11) as desired, we now have 
fT T - T T - T V T f 8 x] 


n r T-T T - T . V_ x 

X ♦ ‘ • 8 y 

" 1 *8 


+ ; {< T x. + *„> C<M “*1 + (T yy + V co * a 2 + T zz c08 °3 ' T vl 
+ n {- F *H co,a l - F AA COBa 2 + slDEco8a 3^ J dt 
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vac 

With the exception of s, each of the other terras in (12) is negative 

m 

on the average and represents a loss. We can therefore define the following: 


ft 


ideal vel 


b T 


vac 


s dt 


m 


(F-13) 


back pressure loss 


T - T 
s vac 


m 


-)dt 


(F-14) 


gimbal loss = 


T - T 

( - y - ~ 8 - ) dt 

m 


(F-i5) 


gravity loss = 


turning loss 5 



tb V 


— 

I . 

g„ 


c a V I 

y 

8„ 

1 

f'b r 

zj 


dt 


1 

m 


(F-16) 


(T + T )co8a, + (T + T )coBa„ + 
xs xx 1 yy ys 2 


T COSO - T 

2Z 3 A v dt (F-17) 


aero loss 


D 

1 

t ' 10 

I 

” F AN COBa I ” F AA CO80l 2 + SIDEco8tt 3 


4 




characteristic vel. “ ideal vel - back pressure loss - gimbal loss 
These losses are such that 


= ideal velocity - pressure loss - gimbal loss - gravity loss - turning 
loss - aero loss. 
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The definition of ideal velocity used here has the added advantage of being 
analytic in the atmosphere. 

However, if other definitions of characteristic velocity are desired, they 
can be calculated easily from (13) through (15), i.e. 


ideal velocity - pressure loss 


ideal velocity - pressure loss - gimbal loss 

* 'a 

It is suggested therefore that (13) through (18) be used in RAGMOP for the 
ideal velocity and losses integration. 


r t. 


T 

-i dt - 
m 


T 

— dr- 
ill 
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Appendix G 

PARTIAL DERIVATIVES OF THE PAYOFF AND CONSTRAINTS WITH RESPECT 

TO THE FINAL STATE* 

The partial derivatives of the various payoff and terminal constraint 
quantities with respect to the final state (position xyz, velocity xyz, mass m) 
are sometimes required to determine the influence coefficient of Paragraph 
3.3.4. These derivative calculations are presented in this appendix. 

Mass M (final cutoff weight) 


8m m 3m 
3 x 3y 


„ o 


3m m ^ 


Inertial Velocity V, 


3V I 3V I 


3V t 

_1 , 

• V 
3y I 


I z 

“ “ v" 


Extracted in large part from reference 1 (see Section IV). 
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|f - A » 63 + B(V s » 13 - W s N 33 ) ♦ |f 


| - A V + B( Vl4 - H sV + 


7 ’ A "65 + B(V s"l5 - V 3 5> + 7? 

3y 3y 
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*7 - * *66 + “< V s N 16 - V»> + ^ 

dz Bz 


if - 0 


Flyback Ranee R, 


3Rp 3R p 3Rp 3 R f 3Rp 3Rp 3R p 


3x 9y 3z 


3x 3y 3z 


The Ug, Vg, and used in the above derivatives are the velocity components 
in the spherical inertial coordinate system. 


Itte used above are defined as fellows: 

N is a 6 x 6 matrix which can be partitioned into four 3x3 submatrices 


K N 

L b iii I "ivj 


The four subraatrices are presented below. 


Defining 


6 1 = 2 - 9 0 


where 


0q - launch site geodetic latitude 


0 ** cos 


, x cosA sinQ, + y cos0, 
-1 z 1 ^ 1 


z sinA 2 sinO^ 


where 


A = launch azimuth 
z 


G-ll 








W s (N 3 i cos6 + N 21 sin6) 
r sin6 

W S^ N 32 COsS + N 22 sin6 ' 
r sine 

W g (N 33 cose + N 23 sine) 
r sine 



’i 


A 
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(V - N >' s > /> 


3x ■ lN n "s “““ - a 3i V /r 


-if ■ (K 12 u ? Ct " 9 - "32 ' J S )/r 


ST - <*13 - s ctne - "33 V /r 



~sin 

A 

2 

cos A 

2 

sine, 

-cos A 

A = 

0 


cosO^ 


sin0^ 


_cos 

A 

z 

-sin A 

2 

sin0^ 

sin A 


z 1 


z 1 


N =0 
III 


11 

r sin6 


12 

r sin0 


: j2 

r sin0 


Thus, N is the matrix of first partial derivatives 

M JS 

w 3P 

where S is the 6x1 vector of spherical coordinate system state components 
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'+» 

r 

^ & J 

and P is the 6x1 vector of plumbline components 




